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Probabilistic Finite Element Methods 


Probabilistic finite element methods (PFEM), synthesizing the power 
of finite element methods with second-moment techniques, are formulated 
for various classes of problems in structural and solid mechanics. 
Time-invariant random materials, geometric properties and loads are 
incorporated in terms of their fundamental statistics viz. second- 
moments. Analogous to the discretization of the displacement field in 
finite element methods, the random fields are also discretized. 
Preserving the conceptual simplicity, the response moments are 
calculated with minimal computations. By incorporating certain 
computational techniques, these methods are shown to be capable of 
handling large systems with many sources of uncertainties. 

By construction, these methods are applicable when the scale of 
randomness is not very large and when the probabilistic density 
functions have decaying tails. The accuracy and efficiency of these 
methods, along with their limitations, are demonstrated by various 
applications. Results obtained are compared with those of Monte Carlo 
simulation and it is shown that good accuracy can be obtained for both 
linear and nonlinear problems. The methods are amenable to implementa- 
tion in deterministic FEM based computer codes. 
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CHAPTER 1 


INTRODUCTION 

Traditionally, engineering analysis has been based on deterministic 
models with well-defined parameters. However, it is increasingly being 
recognized that uncertainties are often associated with parameters such 
as material and geometric properties, forces and boundary conditions and 
that these should be adequately modeled. An example is the degradation 
of material properties with time as a result of fatigue, wear and long- 
term creep; such changes in material properties can be treated as 
uncertainties. In general, the random uncertainties which are included 
in a stochastic process can be classified into three major categories: 
(1) physical uncertainty, (2) statistical uncertainty and (3) 
uncertainty in the model. A detailed discussion of these topics can be 
found in, for example. Refs. [1-4]. Theoretically, these uncertainties 
can be modeled as random variables or random fields governed by joint 
probability density or distribution functions. In practice, the exact 
joint probability density functions are not always available; it is more 
likely that only the first few moments such as the mean and covariance 
are known. 

Uncertainty analysis in structural mechanics has concentrated on 
problems of an almost totally stochastic nature. Within this setting, 
even a single degree of freedom system with nonlinearities poses a 
formidable challenge and has not been solved satisfactorily. The most 
commonly employed solution technique is Monte Carlo simulation (see e.g. 
[3]). In general, these simulation procedures are computationally 
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expensive, even though they are easily applicable to both linear and 
nonlinear systems* For linear systems, nonstatistical methods such as 
"second-moment analysis", are available [2]* A related second order 
perturbation technique applied to a special class of linear structural 
vibrations is discussed in [5]* The emphasis is on the modal decoupling 
of the equations of motion with uncertain damping* The "second-moment 
analysis" has also been extended in [6] to define the mean and variance 
of vector functions* This formulation is mathematically elegant and 
Kronecker algebra and matrix calculus are employed* While this 
formulation has also been extended in [7] to linear stochastic systems 
with colored imiltiplicative noise, the direct application of this 
technique to nonlinear structural dynamics is not feasible, because in 
most nonlinear structural analysis, concern lies more with deviations in 
loads from a deterministic path and in uncertainties in material 
properties to which a value can be assigned rather than completely 
stochastic loads or systems. 

The research reported here can be subdivided into two parts: (1) 

development of a variational principle to embed the probabilistic 
character of the constitutive properties and loads (which are part of 
the boundary conditions and body forces) and to obtain the corresponding 
probabilistic character of the nodal forces; (2) the determination of 
the probabilistic distribution of the response (displacement and stress) 
from the probabilistic description of the nodal forces* 

The main thrust of this research has been to integrate second- 
moment based techniques with finite element methods in a method called 
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probabilistic finite elements (PFEM); finite elements are currently the 
most versatile tools of analysis in large-scale structural and solid 
mechanics. Through this synthesis, we have developed versatile and 
efficient techniques for probabilistic analysis. The investigation is 
restricted to time invariant uncertainties which may be present as 

t 

discrete random variables or random fields in material, geometric 
properties and/or forces. These methods are applicable when the 
uncertainties are not very large and when the probabilistic density 
functions or histograms have decaying tails. The most appealing 
features about PFEM are its conceptual simplicity, ease of computer 
implementation and the flexibility to accommodate efficient numerical 
techniques at every stage of the methodology. 

PFEM has been formulated for various classes of problems in 
structural and solid mechanics. In the next chapter, methods are 
developed for nonlinear structural dynamics with discrete random 
variables. In Chapter 3, random fields are modeled, essentially by 
discretization. The encumberances of correlated random variables are 
avoided by an eigenvalue transformation to the space of uncorrelated 
random variables. In Chapter 4, these methods are derived from 
variational principles. The linear formulation is obtained from the 
potential energy variational principle and the nonlinear counterpart is 
derived from the principle of virtual work with appropriate stress and 
strain measures to account for the large deformation. In Chapter 5, 
numerical applications in elastoplastic mechanics are studied in detail, 
along with improved computational techniques. The summary and 
conclusions are presented in Chapter 6. 



CHAPTER 2 


PROBABILISTIC FINITE ELEMENTS FOR NONLINEAR STRUCTURAL DYNAMICS 
2.1 Introduction 

It is important to be able to treat the effects of uncertainties in 
a reasonably economical manner; standard Monte Carlo procedures are 
simply too expensive. Furthermore, the methods should be designed so 
that they ean be incorporated into widely used finite element programs 
in a natural and concise manner. Thus, the approach should be 
integrable with the elemental discretization and nodal assembly 
procedures that characterize finite element theory and software 
implementation. 

In the next section, the formulation of the probabilistic finite 
element method (PFEM) is presented. The method is applicable in 
structural dynamics with discrete random variables with or without 
correlation. In Section 2.3, the computational aspects of PFEM are 
discussed. In Section 2.4, the analysis of a two degree of freedom 
spring-mass probabilistic system is then given. Results are also 
presented for a ten-bar probabilistic system with nonlinearities. The 
proposed PFEM method is compared to (1) Monte Carlo simulations (MCS) 
and (2) Hermite-Gauss Quadrature (HGQ) schemes. All these methods are 
schematically depicted in Fig. 1, highlighting the major computational 
steps. In Section 2.5, the relative performance of the PFEM as compared 
to the other two methods is discussed. The reason for the limitation of 
each solution technique is also presented. 
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2.2 Formulation of the Probabilistic Finite Element Method (PFEM) 

We consider the structural system to be governed by the following 
system of nonlinear algebraic equations which arises from a finite 
element discretization: 


M d + f(b, d, d) - F(t) (2.2.1) 

/w rv m m /v /w 

where M, f, d and F are the generalized mass, internal force, 

/v /v /v /w 

displacement and external force respectively; and a superscript dot 
represents material time (t) derivative. While the Internal and 
external nodal forces are obtained from one variational statement, they 
are segregated for convenience. The probabilistic effects are described 
through the q-dimensional random vector b; this can Include the 
probabilistic distributions of the material properties; the mass M is 

r*s 

assumed to be deterministic. All these probabilistic distributions, as 
reflected in the variance of the material properties, the composite load 
spectra, etc. are represented by the generalized variance vector, 

Var(b). We shall denote the expected value operator by E[ ] and use 
second order expansions so E[ ] is given by 

- i 

* i jb7Tb- Cov(b i- V (2 - 2 - 2) 

where <j> is a vector function of the random variables. The superposed 
bar denotes "at the mean value of b” and the symbol Cov represents the 
covariance; summations on i and j from 1 to q are assumed. If b^ is 
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uncorrelated to b. 

J 


for i * j, then 


Cov(b^ t bj) ■ 0 for i t j 


and 


CovCb^, b^ - Var(b^) no sum on i 


Applying the expected value operator to Gq. (2.2.1) yields 


'tail * E[«cb, d. i>] ■ 


Employing Eq. (2.2.2) and the chain rules: 


*t«sl ■«5 + i«5T[ib- Cov(1> i' b j ) 


3^ 


<*] ■ i 


♦i 


3 2 f 

a^ab 


Cov(b it bj ) 


3C 3v 3K 3d 

+ t' 3b^ 3bJ + 3b^ 3bJ^ °° v(b l ’ b j J 
2** 

4(£TE^ + £3^fb-l V 


(2.2.3a) 


(2.2.3b) 


(2.2.4) 


(2.2.5a) 


(2.2.5b) 



and 
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E[ F(t )] » F(t) (2.2.5c) 

where cl and d have been replaced by v and a respectively. The C and K 
matrices are the damping and stiffness matrices, respectively. They are 

3f 

C « — (2.2.6a) 

~ 3v 


and 


3f 

K - tt (2.2.6b) 

~ 3d 

In the case of a linear structure, f is given by 

f - C v + K d (2.2.6c) 

A# ^ «V V 

For simplicity, let us assume that Eqs. (2.2.3) holds. This 
assumption is quite suitable for finite element models which are built 
up from discrete structural elements, such as bars and beams. Using 
this simplification and applying perturbation techniques on Eqs. 

(2*2.5), Eq. (2.2.4) can be shown to yield 

Ma+f « F(t) 

(V 


(2.2.7a) 
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and 


M Aa + C Av + K Ad ■ AP 


(2.2.7b) 


where 


3C 3v 3K 3d 


AF 


- l ml Itsbj »bj + 3 bj 3bj ) »“<V 


a 2 ? 


+ T-r Var( V> 

3b j J 


(2.2.8a) 


Aa 


~ 2 


. 1 S 2 ! 

f I -5v»r( b ,) 


j-1 3b: 


j 


(2.2.8b) 


- _ 1 


2 — 
q 3 v 


4 ,._£_ Var(V 


(2.2.8c) 


and 


A d 


— i i 

i'T.I 


2— 

q 3 u 


2 2 Var(b 1 ) 
j-1 3bj J 


(2.2.8d) 


Once Aa, Av and Ad are obtained by solving Eq. (2.2.7b), the second- 
order means are 



Ef al a a + Aa 


(2.2.9a) 



9 


E[ cl] » E[v] a v + Av 


(2.2.9b) 


and 


E[d] - d + Ad (2.2.9c) 

If one is interested in the deviations in response from a 
deterministic path due to the uncertainties in material properties to 
which a value can be assigned, the number of time integrations 
(simulations) reduces to only two. These two simulations are 

M a + 7 « F (2.2. 10a) 

w /v ^ m 


and 


M Aa + C Av + K Ad - AF 

IV /V* IV «V V 


(2.2.10b) 


where 


Aa 


q 3a 

l inr Ab 

j-i 3b j 


j 


Av 


q 3v 

3 Z -l ^ ^ 


(2.2.10c) 


(2.2. lOd) 



Ad 


AF 


q 3d 
q 3f 

-I (,r 

j-i j 


j 


) »», 

d, v » constant J 
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(2.2. lOe) 
(2.2. lOf) 


and Ab is equal to the preassigned value of b . It is also necessary 
J J 

to obtain the sensitivity vectors, i.e. 3( )/3b^; see Section 3. 

Finally, once the means and the sensitivity vectors are determined, 
the variance vectors can be computed easily by the following first order 
formulas 


3a 3a 


v. r (.) - CovO^, bj) 


(2.2.11a) 


3v 3v 


V.i-(v ) . I (^-)(jr-) covCbj, bj) 

l »J*i A J 


(2.2.11b) 


3d 3d 


Var(d) . I (jirKaiB Co ''< b 1 . »j> 

1 J 


(2.2.11c) 


In the case of uncorrelated b’s, the covariance matrix becomes a 
diagonal matrix and the diagonal terms are denoted by 


Var( 


q 3a 


* J-l Var( V 


(2.2.12a) 


Var( 


q 3v _ 

~ -.2 


Z> - l (ab “J Var < b j> 

j-i j J 


(2.2.12b) 



11 



(2.2.12c) 


Similar procedures can also be developed for the probabilistic 
distributions of stresses. However, this direct approach can be very 
expensive if the number of random variables is greater than the number 
of requested probabilistic distributions of stresses. For this 
situation, an alternative approach, termed an adjoint probabilistic 
stress analysis, is developed. This is described in [12]. 


Remark 1 The uncertainties discussed here are described by discrete 
random variables. Physical parameters, such as material properties, are 
often continuous functions in space. When there are uncertainties 
associated with these parameters, we have random fields. The 
probabilistic distributions at any two points can be represented by a 
"correlation function." One way to adapt the above procedures to 
"random fields" is to first do a finite element discretization of the 
correlation function and thus obtain the covariance matrix. Once this 
matrix is obtained the PFEM method as developed here could be used with 
minor modifications. 


Remark 2 To the author’s knowledge, Eqs. (2.2.7) through (2.2.12) 
represent the first consistently derived second moment probabilistic 
finite element method (PFEM) which can readily be adapted to existing 


deterministic finite element computer programs. The second order 
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terms Aa, Av and Ad are computed directly from the second moment mean 
Eq. (2.2.7b). Consequently Eqs. (2.2.9) are second order accurate and 
Eqs. (2.2.11) are first order accurate. 

Remark 3 The complete probability distributions are not available for 
most random variables except perhaps the first two moments. Methods 
such as MCS or HGQ usually require knowledge of probability density 
functions. The PFEM method requires only the first two moments and is 
therefore widely applicable. 



2.3 Computational Aspects of PFEM 
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The computing procedures essentially involve time integrations of 
the various equations derived in the previous section. In general, the 
sensitivity vectors can be obtained directly by integrating the 
sensitivity equations in time. However, this is not possible for some 
nonlinear systems. In such cases, the usual procedure is to calculate 
the derivatives by finite differences [I]. Calculating the finite 
difference derivatives increases the computation for a probabilistic 
system. However, results obtained are excellent when compared to the 
solutions obtained by other methods. The computing procedures for 
linear and nonlinear systems are described separately below. 

Linear Systems 

For a linear system, Eqs. (2.2.7) become 

M a + C v + Kd - F(t) (2.3. la) 

/v /v r\j ix 


M Aa + C Av + K Ad - AF (2.3.1b) 

X f*0 IX /V 0*0 

The solutions of Eq. (2.2.7a) and (2.2.7b) are obtained in sequence so 
that the additional computation due to the latter is minimized. The 
solution algorithms, such as implicit and/or explicit time integration, 
used in Eq. (2.3.1a) can be applied directly to Eq. (2.3.1b) with the 
formulation of only one additional vector function aF. 
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If we examine Eq. (2.2.8a) closely, it can be shown Chat aF can be 

_ _ 3v 3d 

computed element-wise once v, d, and are given. In addition, 
the corresponding variation of the elemental nodal forces can then be 
assembled into a description of the probabilistic distribution of the 
elemental nodal forces for the complete finite element model. 

It can be easily shown that the governing equations for the 
sensitivity vectors are obtained by differentiating Eq. (2.2.1) with 
respect to b^. They are 


3a 


3v 


3d 


M sr- + c TT“ + Krr 

~ 3bj ~ 3bj ~ 3bj 


3F 

3b 


(2.3.2a) 


where 


3F 

3b! 


3f 

3b: 


(2.3.2b) 


d,v 


constant 


or 


3F 

/Ni 

abT 


3C 


3K 


3bj ~ 3b j ~ 


(2.3.2c) 


From Eqs. (2.3.1a-b) and (2.3.2a); it can be seen that the whole 
procedure uses the same effective stiffness matrix so only one matrix 
needs to be triangulated. 

To evaluate the mean and variance from Eqs. (2.2.9) and (2.2.12), 
the total number of time integrations required is q + 2. These are: 
one integration to evaluate the displacement, velocity and acceleration 
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at the mean value of b (Eq. 2.3* 10a); f q f integrations to evaluate the 
sensitivity vectors (Eq. 3.2a); and one more integration to evaluate the 
second order variations (Eq. 2.3.1b). The computational steps involved 
in PFEM are shown in Fig. 2. Notice that all time integrations employ 
the same effective stiffness matrix; parallel computation procedures 
could be employed, thereby increasing the efficiency tremendously. 


Systems with Material and Geometrical Nonlinearities 

As in the linear case, the displacement, velocity and acceleration, 
at the mean value of b is obtained by integrating Eq. (2.2.7a). The 
relative merits between implicit and explicit time integrations are 
considered here for a probabilistic nonlinear system. 

By total differentiation of Eq. (2.2.10a) with respect to Bj , i.e. 
d/dbj , we have: 




(2.3.3a) 


and 


d 2 i 


d 2 T 


M — + — 

~d b j d bj 2 


Equations (2.2.7a) and (2.3.3a) can be written as 


(2.3.3b) 


~ ~n+l + ~n+l " ^n+l^ 


(2.3.4a) 
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and 


M ^±L + K 3 ~ n+1 


~ 3b. 


3b, 


3f 

“ ibj' ~n+P 


(2.3.4b) 


where £ n+ ^ and K are the internal force vector and the "tangent 

stiffness matrix", respectively, evaluated at T>, d" , and t ,. 

— -n+1 n+1 

Equations (2.3.4a) and (2.3.4b) can be solved by the implicit Newmark-g 
algorithm [10]. The "mean value" equation (2.3.4a) can be solved by 
Newton-Raphson iteration 


r* *fl v+ 1 * T v 
— —n+1 — a+1 


(2.3.5a) 


where the residual vector Is given by 


-n+1 **— n+1 —n+1 


[F - f , - M a ,1 V 
L ~n+1 ~n+l ~ ~n+l J 


(2.3.5b) 


and the effective stiffness matrix is 


2 7|V 


K * {M + gAt Kl 

/v I*** 


(2.3.5c) 


The symbol v represents the equilibrium iteration counter at time step 

n+1 and iterations are repeated until Aa v *| approaches zero. 

-n+l 

Similarly, the first order sensitivity equation (2.3.4b) can be 
written as 



a a-i 

3b, 
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(2.3.6a) 


where 


3d . 3d 3v 3 a 


(2.3.6b) 


It is observed here that the effective stiffness matrix K is identical 


in both Gq. (2.3.5a) and Gq. (2.3.6a). Since the triangulated K is 

3 ?n+l 

given during the iteration procedures, can be obtained simply by 

3 j 

forward reductions and back substitutions; therefore, the number of time 
integrations is still q + 2. 

The main advantage of employing implicit time integration is its 
unconditional stability. Therefore, the above methods are best suited 
for structural dynamics problems dominated by low frequency response. 

For Impulsive and short duration transient problems, Gq. (2.2.7a), 
(2.3.3a) and (2.3.3b) can alternatively be solved by explicit 
integrations. Since f(b, d) is nonlinear, the sensitivity vectors can 
be obtained by central-differences. Equations (2.3.3a) and (2.3.3b) are 
approximated by 


+ 

a - a 


f + - f- 


S (=SSjH - (=S5f-) “ o 


(2.3.7a) 


and 



2a + a 
** ** 



') + 



- 2f + f 



0 

** 


(2.3.7b) 
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where 


a ■ a(b + Ab, ) 

/V /V /V 1 

a =* a(b - Ab, ) 

l*W /V IV ^1 

and 


f + - f (b + Ab. , d + ) 

/V V /X IVJ V 

f" - f(b - Ab, , <f) 

<*w IV IV ivl V 


d and d are similarly defined and Ab, is defined by 

#v 


Ab , * (0 > Oj • • • t Ab . ) 0, • • • 9 0) 

~j j 


(2.3.7c) 

(2.3.7d) 

(2.3.7e) 
(2.3.7f ) 

(2.3.7g) 


where T denotes the transpose. With this computational procedure, the 
total number of time integrations would still be q + 2. However, the 
number of Internal force calculations would be 2q + 1. These are: one 

integration for the mean Eq. (2.2.7a) and 2q integrations with finite 
differencing for Eq. (2.3.7a) and Eq. (2.3.7b). Apart from purely 
implicit or purely explicit algorithms, mixed time implicit-explicit 
algorithms [4] could also be employed so that the attributes of each of 
the algorithms can be achieved. 
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2.4 Numerical Examples 

Example 1 ; A Two-Degree-of -Freedom Spring-Mass System. 

The performance of PFEM, the new method developed here. Is 
evaluated via a two-degree— of -freedom spring-mass system. The mean is 
second order accurate and the variance is first order accurate in this 
example. The computed results are compared with those obtained 
employing (1) Monte Carlo Simulation (MCS) and (2) Hermite-Gauss 
Quadrature (HGQ) schemes. The two latter methods as implemented here 
are reviewed in Appendix A. 

The problem statement is depicted in Fig. 3. A sinusoidal vector 
forcing function is used: 


F(t) 


0.0 


25.0 x 10 6 sin 2000t 


(2.4.1) 


The random spring constants Kj and K 2 are normally distributed with a 
coefficient of variation (i«e« a/y) equal to 0*05« The mean spring 
constants are 24 x 10 6 and 12 x 10 6 respectively. The deterministic 
masses and m 2 are 0.372 and 0.248 respectively. A stiffness- 

proportional damping of 3 Z is included. The probabilistic equations 
derived earlier are solved by the implicit Newmark-g method [3]. The 
mean amplitude d^ is depicted in Fig. 4, for all the three numerical 
methods — PFEM, HGQ and MCS. The PFEM solution compares very well with 
the other two methods. For the variance of d^ the PFEM solution plotted 
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in Fig* 5, seems to overshoot the variance at large times* The mean and 
variance of are similarly compared and depicted in Figs. 6 and 7. 

The maximum coefficient of variation of the displacements d^ and d^ are 
found to 0*13 and 0*10 respectively* The ±3o bounds for the 
displacements d^ and d^ are plotted in Figs* 8 and 9 respectively* 


Example 2 s A Ten-Bar Probabilistic System with Material and Geometrical 
Nonlinearities * 

The problem statement is depicted in Fig* 10. The load time 
function, which is also shown in Fig* 10, is applied at node 3. This 
particular load-time history is chosen such that only four of the ten 
bars, elements 1, 3, 7 and 8, will yield. Therefore the probabilistic 
model can be simplified by choosing the yield stresses of these four 
elements as the normal random variables which have the major impact on 
the response* The coefficient of variation is 0*05* Since the other 
six elements do not come close to yield, they are considered 
deterministic variables* With this approach, instead of 59049 analyses, 
only 81 analyses are required for the Hermite Gauss Quadrature method. 
The justification for this drastic simplification is explained in detail 
in Appendix A. 

For the PFEM method, the finite difference derivatives are 


evaluated with an interval equal to 0*05 b^ and the equations are 
solved by explicit time integration* The mean is second-order accurate 


whereas the variance is first-order accurate* The Monte Carlo 


Simulation results are obtained with 400 simulations. 


The probabilistic displacement and stress solutions at selected 
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locations are given in Figs, 11 through 14, The maximum coefficient of 
variation of the displacement of node 1 is found to be 0.13 and that of 
the stress in element 1 is 0.11. For this example, the three methods 
(PFEM, HGQ and MCS) have been employed and they all compare quite well. 

The bounds of the displacement and stress can be estimated based on 
the Chebyschev inequality 

P( | x - y | > no) < , n > 0 (2.4.2) 

n 

2 

where y ■ E(x) and a ■ Var(x). The ±3o bounds (i.e., n ■ 3) for the 
displacement and stress are plotted in Figs. 15 and 16, and the 
solutions can be expected to be within these bounds with 89% confidence 


level 
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2.5 Comparisons Among the Three Methods and Conclusions 

Based on these numerical studies , we have drawn the following 
tentative conclusions: 

1) Although all three methods agree very well and are evidently 
comparable in accuracy, PFEM is the most efficient solution procedure 
for small to medium size problems. The relative computational 
efficiency of the three methods is summarized in Figure 17. 

Relative Computational Efficiency of Three Probabilistic Methods. 

PFEM HGQ MCS 

1 8 400 

1 4 60 

Figure 17 

The number of time integrations required for a general structure 
with q random variables can be summarized as follows: 

i) PFEM with partial derivatives evaluated directly: q + 2 

ii) PFEM with partial derivatives evaluated by finite difference: 2q 
+ 1 

q 

iii) HGQ with three-point quadrature : 3 

iv) MCS with simple Monte Carlo Simulation of sample size N : N 


2 bar 
Structure 

10 bar 
Structure 
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2) Although PFEM is expected to be most accurate when the variances are 
small, it performs quite well even when the response shows a large 
coefficient of variation (e.g., 0.13 for the displacement at Node 1 in 
the ten-bar structure). This could be attributed partly to the nature 
of the probabilistic distribution. For most distributions, values of 
response far- away from the mean are less likely to be found than those 
near the mean. Hence second moment analysis about the mean turns out to 
be quite accurate. 

3) The three methods are applicable to linear and nonlinear systems. 

In linear systems the partial derivatives can be obtained directly. In 
nonlinear system the brute force method is to obtain these derivatives 
by finite differences. We are currently investigating ways to compute 
these derivatives efficiently. However, the methods are problem 
dependent. 

4) A minor drawback of PFEM is that its accuracy deteriorates for large 
times even with structural damping. An explanation is given in Appendix 
B. We are currently investigating several ways of improving this. 

5) PFEM can be easily incorporated into widely used finite element 
programs. 

6) A PFEM analysis can be obtained with q + 2 simulations if CovCb^ 
b j ) * 0 for i * j. For this purpose the b^ must be transformed into 
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anocher set of random variables through an eigenproblem such 
chat Cov(c^, c^) * 0 for i * j and in most cases only a few modes are 
sufficient [13]. 


7) Currently the PFEM is being extended to the transient analysis of 
nonlinear continua. The details of the method can be found in [5]. 

Since this method involves only matrix and vector assembly it can 
be incorporated in a natural and concise manner in general purpose 
finite element programs. 
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Schematic of Probabilistic Methods 
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A Two Degree of Freedom Example 
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CHAPTER 3 


RANDOM FIELD FINITE ELEMENTS 

3.1 Introduction 

At the present time, probabilistic methods in mechanics, for 
problems involving time-independent uncertainties, can be broadly 
classified into two major categories: (1) methods using a statistical 

approach and (2) methods using a non-statistical approach. The 
literature in these areas is quite considerable and so only a few sample 
references are Indicated below. 

Simulation, involving sampling and estimation, is the most 
prevalent statistical approach. Direct Monte Carlo simulation, 
stratified sampling and Latin Hypercube sampling are some of the 
frequently employed simulation techniques. A comparative discussion of 
these techniques can be found in, for example. Refs. [1,2, 3, 6]. These 
techniques, however, have their limitations. Transformations of the 
distributions are necessary before simulation can be done [4, 5, 8, 9]. 

This implies, of course, that the multivariate distribution function 
needs to be known for simulation. The topic of transformation 
techniques is still an area of current research. Furthermore, since the 
accuracy of sampling techniques depends on the sample size, in 
accordance with the "Weak Law of Large Numbers" [4,11], simulations can 
become prohibitively expensive; hence the interest in non-statistical 
methods. 

Non-statistical approaches include numerical integration 
[10,15,16], second-moment analysis [4,5,7,9,14-17], and stochastic 
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finite element methods [11-13,15,16]. Particularly, second-moment 
techniques have proven to be accurate and efficient in structural 
mechanics. A major advantage of these techniques is that the 
multivariate distribution function need not be known but only the first 
two moments. An inherent limitation of second moment analysis is that 
the uncertainties cannot be too large, i.e., variances of the random 
variables cannot be large when compared with their mean values. 
Typically, the maximum coefficient of variation is around 10% although 
it has been shown that it could be as high as 20% for acceptable results 
to be obtained [5,14]. 

Linear problems in structural mechanics with uncertain parameters 
have been solved by second-moment analysis [11-13]. However, similar 
solution techniques for nonlinear problems in structural dynamics are, 
to the authors' knowledge, nonexistent. Recently, the authors have 
developed probabilistic finite element methods for nonlinear structural 
dynamics [15,16]. The methods are applicable to correlated and 
uncorrelated discrete random variables, though they are limited to 
discrete structures such as spring-mass systems and nonlinear truss 
structures. 

The herein proposed method is applicable to nonlinear structural 
dynamics problems with random fields — both homogeneous and 
inhomogeneous. In the next section, the formulation of the 
probabilistic finite element methods for linear continua is outlined. 

In Section 3.3, the procedures for the transformation of the full 
covariance matrix to a diagonal matrix are discussed. In Section 3.4, 
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the computational procedures using the transformed random variables are 
given. The PFEM, as applied to continua with material and geometrical 
nonlinearities, is formulated in Section 3.5. Applications to a one- 
dimensional elastic/plastic wave propagation problem and a two- 
dimensional plane-stress beam bending problem are described in Section 
3.6* The results and conclusions are presented in Section 3*7. 
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3.2 Probabilistic Finite Element Methods (PFEM) for Linear Continua 
The linear finite element equations ares 


K d - F (3.2.1) 

/V SSJ 

where the stiffness matrix is 

K - f B T D B dQ (3.2.2) 

Q A# IV 

The transpose is designated by a superscript "T"; the generalized 
gradient matrix, material response matrix, nodal displacement vector and 
nodal force vector are denoted by B(x), D(x,b), d(b) and F(b), 

/v /w <v /w /v 

respectively; x are the spatial coordinates; ft is the domain and b(x) is 
a random function. In this formulation, b(x) can be a random material 
property or a random load. 

The basic idea of the "Second Moment Analysis" in PFEM is to 
expand, via Taylor Series, the d, D and F matrices about the mean value 

of b and to retain only up to second order terms. Equations will then 

be obtained for the mean values of the nodal displacements and the 
covariances of the nodal displacements in terms of the derivatives of 
the nodal displacements with respect to the random variables. 

Similarly, the mean and covariance of the element stresses and strains 
are obtained. 

The random function b(x) is approximated using shape functions 


^(x) by 



q 

b(x) - Z N. (x)b 
~ i-1 1 ~ 1 

where b, are the nodal values of b(x), that is the values of b at x. , 

A ^ 'vl 

i * 1 1 • • • i • 

To derive the PFEM matrix equations, the following notation will be 
used. For a given function g(b) and a snail parameter e: 

b(x) ■ E[b(x)J mean value of b, i.e. the 

expectation E[ ] of b(x) 

db^ » eAb^ « eCb^ - b^) first order variation of b^ 

about b^ 

2 

dbjdb^ * e AbjAb 2 second order variation of b^ and 

b 2 about b^ and b 2> respectively 
g(x) * g(x,b(x)) value of g evaluated at b 

partial derivative of g with 
respect to b^ evaluated at b 
mixed partial derivatives of g 
with respect to b^ and b£ 
evaluated at b* 

The random function is defined by its expectation b(x), coefficient 
of variation a and autocorrelation R(b(x^) ,b(xj)). The mean and 
variance are approximated by the same shape functions as b, so 


7 -il_ 

*1 3b l 


g 


g 


3 2 g 

b l b 2 3b l 9b 2 


45 

(3.2.3a) 


q 

E[b(x) ] =» Z N (x)E[b ] 
~ 1 ~ x 


q 

Var(b(x) ) = Z N.(x)N (x) Cov(b ,b.) 

< 4_1 1 ~ J ~ 1 J 


(3.2.3b) 

(3.2.3c) 



where 
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about b via Taylor series: 

3 . ( 2 . 4 ) 

( 3 . 2 . 5 ) 

( 3 . 2 . 6 ) 

.2.6) into Gq. (3.2.1) and equating 
and second order equations 

( 3 . 2 . 7 ) 

( 3 . 2 . 8 ) 



i * 1 f «««} q 


( 3 . 2 . 9 ) 
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(3.2.10) 



F. - / B T D, Bddii 

D. ' ~ ~D ~ ~ 

i 0 i 


Second Order 



terms) 



(3.2.11) 


where 


- 1 q - 

d, - -x I d, . Ab.Ab 

~ 2 2 t,j-i ~ b i b j 1 1 


(3.2.12) 


and 


h- 2 - il/ b 2 i <«»] 

i.j-1 i j n i j 


- [/ l T S bl » 3 b , dc ]l'‘ b 1 4 bj 


(3.2.13) 


Once d, d, and 7. are obtained by solving Eqs. (3.2.7), (3.2.9) and 
(3.2.11), respectively, the mean and autocovariance matrices for the 
nodal displacement are given by 


00 

E[d] * / d(b) f(b) db 

/W » (V A# Af /W 


(3.2.14) 


and 


Cov(d i ,d^) =■ f (d 1 - d 1 )(d j - d j )f (b)db 

AW • A> M 


(3.2.15) 
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respectively, where f is the joint probability density function, d* is 
the ith degree of freedom of d, and b is the random-variable vector 

a# 


b 


(b ^ • * • > 



(3.2*16) 


The second order estimate of the mean value of d is obtained by 
employing Eq. (3.2.5) in Eq. (3.2.14) to give 


E[d] 

<v 



ib 1 b j Cov(b i- b j ) t 


(3.2.17) 


The covariance, Cov(b^ ,bj ) is obtained from the given expectation 
E[b(x)], coefficient of variation a and autocorrelation R(b(x^) ,b(xj)) 
as follows 


Cov(b^,bj) * (Var(b(x i ))Var(b(Xj))] 1(/2 R(b(x i ),b(Xj)) (3.2.18) 


where 


Var(b(x t )) - a 2 E[b(x i )] 2 


(3.2.19) 


Similarly, the first order accurate Cov(d*,d^), which is consistent with 
a second moment analysis, can be shown to be 


Cov(d i ,d^ ) - E df d? Cov(b ,b ) 
, D D r s 

r,s-l r s 


(3.2.20) 
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The strain and stress vectors for a typical element "e” are 


e - B(x) d® 

AW aw aw 


(3.2.21) 


and 


o - D(x,b) e 

/W A# AW 


(3.2.22) 


0 

where d is the element nodal displacement vector. Since B is a 
deterministic function of x, the mean value and autocovariance of e can 

*** AW 

be similarly defined according to Eqs. (3.2.14) and (3.2.15), 
respectively. They are 


E[e] 

AW 



♦ 


{ 

ij- 


5Sib j Cov(b i» b j ) } 


(3.2.23) 


and 


■ 1 E (J'lj ><J f 4b )Tc °'' (b 1 - b J 

i j 


)} 


(3.2.24) 


Employing Eq. (3.2.4) and the element counterpart of Eq. (3.2.5), 
the mean value and autocovariance of a can be shown to be 

AW 

E[o] - D E[e] 


+ ( e [»„ »S. + iSb.b.J ?] Cov(b i- b j , > 

l»j“l i J 1 J 


(3.2.25) 
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and 


Cov(o ,o f ) - { I [(0*8*5? XljVd* ) T 

i»j-l i i 

+ (D? B*5*)(D? B^) 1 + (0*8*5? )(d5 bW 

r* Q^<v /v /v <v <v a# <v »vQ^ /v0^ ** 

+ (D? B*5 e )(D f B f df ) T ]Cov(b. ,b. )} (3.2.26) 

/V 0 , /v M /V »v /vQ J 1 J J 


respectively 
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3.3 Transformation of the Full Covariance Matrix to a Diagonal Variance 
Matrix 

It can be observed from Eqs. (3.2.7) through (3.2.13) that the 

determination of d, <T. and T, involve one factorization of K and q + 2 

forward-reductions and back-substitutions. The latter operations 

consist of one solution to evaluate d (Eq. (3.2.7)); q solutions to 

evaluate T. (Eq. (3.2.9)) and one more solution to evaluate d„ (Eq. 

~b . 

i 

(3.2.11)) where d„ la obtained by multiplying the joint probability 
density function with Eq. (3.2.11) and integrating over the domain 
of b to yield 

Kd'-?’ (3.3.1) 


where 


V 



and 



^b^w 


(3.3.2) 


ii- ; (iSb.b. 

i»j-l i j a i j 


U J T J b 5 3b da ]t Cov<b i- b J ) 

a i j J 


Hence, from Eq. (3*2.17), the mean value of d is simply 


(3.3.3) 


E(d] - d + d' 


(3.3.4) 
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Even though the above computations are compatible with the elemental 
discretization and nodal assembly procedures that characterize finite 
element theory and software, the number of matrix multiplications is 
proportional to q(q + l)/2. This would be unacceptably expensive. The 
large number of computations arise from the double summations in 1 and j 
in Eq. (3.3.2) and (3.3.3). To remedy this situation, the covariance 
matrix Cov(b^,bj) is transformed to a diagonal variance matrix 
Var(c i ,Cj) such that 

Var(c i ,c j ) - 0 for i + j (3.3.5) 


and 


Var(c lf cj) - Var(c^) 


for i * j 


(3.3.6) 


Therefore, the number of matrix multiplications is proportional to q. 
The above involves the solution of the eigenproblem: 


G ip » <|> A 

fv rw 


(3.3.7) 


where the G and A matrices are used to denote Cov(b^,bj) and Var(c^,Cj), 
respectively; and ip is a constant q x q fundamental matrix with the 
following properties: 


T T 

" i 


(3.3.8) 
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(3.3.9) 


A “ t S £ 


and 


b - il> c or c * ili^b (3.3.10) 

M A* r*t <-w IV 

1 is the q x q identity matrix and c is a transformed q x 1 random 

V IV 

variables vector. 

With Gqs. (3.3.9) and (3.3.10), the mixed derivatives appearing in 
Eqs. (3.2.21) through (3.2.28), (3.3.2) and (3.3.3) reduce to second 
derivatives. For example, Eqs. (3.3.2) and (3.3.3) become: 

t 2 m \ = ?c c Var(c i ) (3.3.11) 

1 1 i-1 c i c i 1 


and 


i 



1 llJc c - it/ £% c 5 i «] 

i-1 i i a i i 


- [/ 5 T £ Cl 5 ic t d °]J v “< c ±> 


(3.3.12) 


Analogous to modal analysis in structural dynamics problems , only a 
few modes (i.e. Var(c i )) are required to capture the major 
characteristics of the probabilistic distributions. However, the 
highest eigenvalues have to be employed. This is in contrast to the 
modal structural problem wherein the lowest eigenvalues are used. 
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3.4 Computational Procedures Using the Transformed Random Variables c 
Assume that n highest Var(c^) are adequate to characterize the 
probabilistic distribution. The discrete c^ are transformed according 
to 


c i ' Vj 


i a 1 f . . . i n 


(3.4.1) 


and the mean and variance of c are 


E[c] - * E[b] 

M <V /V 


(3.4.2) 


and 


Var(c) * diagonal terms of A 


(3.4.3) 


where 


T 

A * 4 G * 

At 


The zeroth-order matrix equations to be solved are: 


K d - F 

»w (W /v 


(3.4.4) 


(3.4.5) 


The n first-order matrix equations to be solved are: 
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K d 


C 


i 



i - 1, n (3.4.6) 


where 


F _ - F - [/ B T D B d dft] 

~i+2 ~c L J ~ ~c ~ ~ J 

in i 


i * 1 ) . * • i n 


(3.4.7) 


The second-order matrix equations to be solved are: 


K d - F 


(3.4.8) 


where 


f! - I - -kj B^D B d da] 

1 2 . 2 L ^ ~ -c i c r ~ J 


- [/ s\ 5 i c ■iallVarCc ) 
a i i 


(3.4.9) 


It is also interesting to note that Eqs. (3.4.5) through (3.4.9) 
can be put into the general form of 


K d » F 


(3.4.10) 


where K is an (n + 2) x (n + 2) block lower triangular matrix of the 


form 
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K 


K 


K 

: c ] 


0 

K 

0 


0 

0 

K 


• • • 


• • 


K 

~c 


0 0 


K K K 
~ I ~2 


0 

0 

<**# 

0 


0 

/w 

0 

/v 

0 


• • 


• • 


K 0 


K K 


(3.4.11) 


where K, K , K, and K 00 are 
** ^zz 


k - / b t d b da 

<V * /V IV /w 

8 


(3.4.12) 


K - / B T D B dQ 
18 1 


i - 1, ..., n (3.4.13) 


£i - I flS c Vlr < c i>]5 <“> 

8 1 


1. * 1 y • • • | tl 

no sum 


(3.4.14) 


and 


^ I _ u 

K„ J B X [ E D Var(c,)]B d8 
~ 22 2 a 1-1 ~ c i c i 1 


(3.4.15) 


Accordingly, d and F are n + 2 block vectors which are defined by 


d - (d, d , d , d , d ) T 

^ ^ . /vC. ^z 

12 n 


(3.4.16) 


and 
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F 


(F, F , 





( 3 . 4 . 17 ) 


where 


h 


l Q - 

3 Z lc c 

1 i-1 c i c i 


Var(c ± )] 


( 3 . 4 . 18 ) 


The mean and autocovariance matrices for the displacement are: 


E[d] - d + d. 


( 3 . 4 . 19 ) 


and 


CovCd^d 1 ) - S d 1 T L Var(c.) 

j-1 "j C j J 


( 3 . 4 . 20 ) 


The mean and autocovariance matrices for the element strains are: 


E[e] * B a + B dl 

m* Ai rv ^ JL 


( 3 . 4 . 21 ) 


and 


Cov(e_ ,6c) 



( 3 . 4 . 22 ) 


Similarly, the mean and autocovariance matrices for the element stresses 
are: 


(3.4.23) 


+ { I [D Bd* +-J-D B 7*lVar(c,)} 
~c^~ ~c^ 2 ~c^c^~ «* •* i ‘ 


and 


-< 2 .. Je > ■ ( ivmjm? 

+ (D* B e d e )(D f B^d^ ) T 

A# ' ' ^ ^ ' 

+ (D e B e d e VU* B^) 1 

+ c5 i sDa'i^ t ) t ]»«‘i»j (3 - 4 - 24) 

Remark 4.1 While the presentation of the PFEM solution algorithm via 

Eqs. (3.4.10) through (3.4.18) is quite elegant, these equations are not 

employed in practical computations since the formulations of K , K. , 

~c^ ~i 

K 00 and the triangulation of K are unnecessary. Instead, the 
~*22 ~ 

"sequential algorithms" given by Eqs. (3.4.5) through (3.4.9) are 

employed. It should be noted that only the K matrix needs to be formed 

and triangulated and f\ ¥. . „ and 1F„ are obtained by vector computations. 

~ ~i+2 ~2 1 1 " " " 

Remark 4.2 Once 7 is given by Eq. (3.4.5), the d as defined by Eqs. 
(3.4.6) and (3.4.7) are best obtained with parallel computations. 
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Remark 4.3 The matrices D, and D. . can be comoated bv direct 

^ 1 J 

differentiation or by the least-squares fit of the series Eq. (3.2.4) 

around b. Once these matrices are obtained, the transformation Eqs. 

(3.4.1) and (3.4.2) can be employed to yield I) and 1) . Similarly 

_ _ ~ c i ~ c i c i 

F and F can be obtained. 

~ C i ~ c i c i 


Remark 4.4 The PFEM, as developed here, can Incorporate smooth (C°) 
shape functions in Eq. (3.2.4). However, this may result in more 


integration points in the evaluation of K , K. and K„„ in Eqs. (3.4.13) 

~22 

through (3«4«15)« To minimize computations, a super-element that spans 
several elements used for the displacement approximations and in which 
the shape function is a constant, is employed in the numerical examples 


studied here 
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3.5 PFEM for Transient Analysis of Nonlinear Continua 

The transient equations for the finite element model, which account 
for both geometrical and material nonlinearities, are: 


M a(b,t) + f(d,v,b) 

m /w /w 


F(b,t) 


(3.5.1) 


where M is the deterministic mass matrix; f is the internal force 

/v *v 

vector; d, v and a are the displacement, velocity and acceleration 

vectors, respectively; F is the external force; b is the discretized 

<■»* 

random vector and t is the time. Following the same procedures outlined 
in the previous sections, the a, F and f vectors are expanded via Taylor 
series, however, total derivatives are applied to f. The second-order 
formulas for a, F and f are 


a =» a + L a, db . + -r l a db db 

~ ~ i-i ~ b i 1 2 i,j-i ~ b i b j 1 j 


q _ i q _ 

F ■ F + E F, db . + 4- Z F. , db.db. 

~ ~ i-1 ~ b i 1 2 i,j-l ~ b i b j 1 j 


(3.5.2) 

(3.5.3) 


and 


_ q _ 

f » f + z [f. + C v. + K d. ldb 
~ ~ i=l ~ l ~ ~ ~o^ J 1 


+ i,j-l^ 2 ~Vj + 2 ^ ~ b i b j + 2 ~ ~Vj 


+ ^b^ + \^^ db i db i 


(3.5.4) 
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where Che tangential damping matrix and the tangential stiffness matrix 
are defined by 



and 


3f 

sw 

K - rr 
~ 3d 


(3.5.6) 


respectively. Using the above approximations in Eq* (5.1), the q + 2 
solutions for d, v and a are: 

A* 

Zeroth-Order Equation 


M a + f - F 


(3.5.7) 


First-Order Equations 


Ma. + C v + K d - F 


i - 1, ...» q 


(3.5.8) 


and 





i 


i*l» •••> 9 


(3.5.9) 



Second-Order Equations 


M a. + C v- + K d. 



(3.5.10) 


where 


a„ • — 


, <1 _ 

T E K Cov(b ,b ) 

2 i,j-i b i b j 1 J 


, <1 

~o 1 v k k Cov(b. ,b ) 
2 i,j-i ~ b i b j 1 j 


*2 


1 *1 _ 

T 2 dv s Cov(b. ,b, ) 

2 i.j-i % i b j 1 J 


and 


(3.5.11) 

(3.5.12) 

(3.5.13) 


~2 ™ 


q 

£ 


i t2 Vj 


- 2 £b 1 b J 


- v> J ,Cov<b i' b j ) 


(3.5.14) 


The computational effort in solving Eqs. (3.5.8) through (3.5.14) 
can be reduced significantly by transforming the full covariance matrix, 
Cov(b£,bj), to a diagonal variance matrix, Var(c^). Since Eqs. (3.5.8) 
through (3.5.14) are linearized equations, the transformation procedures 
are parallel to those outlined in Sections 3 and 4. If n (recall n < q) 
highest Var(c^) are used, the q + 2 block system becomes an n + 2 block 
system. These n + 2 blocks are: 


i » 1 


» • • • > 


n 
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(3.5.16) 


Ma + C v + K d 

A/ A« A#(2^ A< A#^^ 



M a. + C V 0 + K d. 

ry a« a#^ <v 



(3.5.17) 


where 





i - 1, ...» n (3.5.18) 


and 



1-J 2 ~ C i C i 




V 

~c 


1 


- 5c 1 Sc 1 l Vac(c i> 


(3.5.19) 


Equations (3.5.15) through (3.5.19) are similar to those developed for 
the probabilistic dynamic response of a truss structure with 
uncorrelated random variables [15,16]. Therefore, the numerical 
solution algorithms given in [15,16] can be employed directly for the 
solutions of the above equations. Once a, v, d, a , v , d , a , v 0 

^ ^ ^ a#C^ '**C ^ ^4 ^4 

and d 0 are obtained, the mean and autocovariance matrices can be 
computed according to: 


E[d] - d + d, 

AW A/ 


(3.5.20) 


E[v] - v + v, 

A# »/ *** Z. 


(3.5.21) 


E[a] - a + a 

A# Ay 


(3.5.22) 
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and 


Cov(d i ,d^ ) - { 2 ^ Var(c )} 

r*l r c r 


Cov(v i ,v^) - { l Var(c r )} 

r-1 r r 


Cov(a i ,a^) - { Z a^ a^ Var(c f )} 
r«l c r c r 


( 3 . 5 . 23 ) 

( 3 . 5 . 24 ) 

( 3 . 5 . 25 ) 



65 


3.6 Applications 

The PFEM formulation developed in Sections 2 to 5 has been tested 
by studying two different applications. These are: (1) wave 

propagation in a one-dimensional elastic/plastic bar; and (2) static 
response of a two-dimensional plane stress elastic/plastic cantilever. 

In these numerical examples, the expectation, the spatial 
autocorrelation and the coefficient of variation of the random field 
b(x) are assumed as follows: 

E[b(x.)] ■ b (1.0 + 0x,/L) (3.6.1) 

1 O 1 


and 


R(b(x i > ,b(Xj)) - exp(-|xj-x ^ | /\) (3.6.2) 

ct m 0.10 (3.6.3) 


where x^, L, X, and b(x^) denote the location, the length of the 
bar/beam, the correlation length and the random function at x^, 
respectively. 8 and b Q are constants. It is to be noted that the 
autocorrelation between any two points depends only on the interval 
between these points and not on their locations. The material in the 
bar/beam is assumed to be linear elastic and isotropic hardening, with 
the uniaxial yield stress as a Gaussian random field in the axial 
direction. As can be seen from Eqs. (3.6.1) and (3.6.2), the yield 
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stress is a linear function for the mean and an exponential function for 
the spatial autocorrelation. 

The problem statement for the bar is depicted in Fig. 1. The 
random field is discretized so that q - NUMEL * 32. The probabilistic 
equations derived earlier are solved by the explicit predictor algorithm 
[15] with a slight numerical damping (y ■ 0.55). A near-critical time 
step (At - .000455) is used to keep the number of time steps minimal, 
subject to the stability conditions. 

In the case of the beam, the static response is calculated as a 
function of steadily increasing loading, by an implicit algorithm. The 
random field is discretized with 64 4-node 2D elements so that q * 16 
(NUMEL - 64). 

The mean and variance of the displacement at the free end of the 
bar, the variance and autocorrelation of the stress along the beam are 
computed using PFEM. These results are compared with Monte Carlo 
simulation (MCS) of 400 realizations with a first-order filter [5,9]. 
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3.7 Results and Conclusions 

The mean and the variance of the displacement at the free end of 
the bar, computed by PFEM and MCS, are compared in Figs. 2 and 3. The 
coefficient of variation of the displacement at the free end of the bar 
is found to be ~0.13 (the coefficient of variation of the yield stress 
is 0.10). The PFEM solutions compare very well with the MCS solutions. 
Although both methods compare very well in accuracy, the PFEM in this 
case needs much less computer time than the MCS. The convergence of the 
variance of displacement at the free end of the bar is plotted in Fig. 

4, against the number of modes used in the PFEM computations. It is 
observed that only the largest 8 of the 32 eigenvalues (which correspond 
to the variance of the uncorrelated variables c t ) are sufficient, for an 
error less than 1%. 

The variance of the stress at the fixed end of the beam, with 
increasing loading, is plotted in Fig. 5. The maximum coefficient of 
variation is found to be 0.09, and it occurs when the beam begins to 
yield, at the mean configuration. The PFEM variances are in excellent 
agreement with those of MCS. The autocorrelation of the stress along 
the length of the beam, with respect to the stress at the fixed end, is 
plotted in Fig. 6. As expected, the autocorrelation near the fixed end 
is ~1.0 and beyond that it decreases rapidly. The PFEM autocorrelation 
is in fairly good agreement with that of MCS. As in the case of the 
bar, only a few eigenvalues are found to be necessary. The largest 4 of 
the 16 eigenvalues are sufficient to ensure an error less than 5%. 

Since the random field is handled by discretization, it is easy to 
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Incorporate any specified mean, variance and autocorrelation structure 
in PFEM. As the stiffness matrix corresponding to the mean value of the 
random field appears in all the PFEM equations, the triangulation needs 
to be done only once and the computations are thereby reduced. The 
transformation of the correlated variables to a set of uncorrelated 
variables further reduces the computations as the covariance matrix is 
reduced from a full matrix to a diagonal matrix. However, to do this an 
eigenvalue problem of the covariance matrix needs to be solved. 

Numerical results obtained here suggest that a reduced set of the 
uncorrelated variables is sufficient to predict the response moments 
accurately. The PFEM essentially involve solution of a set of 
deterministic problems, and therefore, they are easily integrable into 
any FEM based code. 
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CHAPTER 4 


PROBABILISTIC FINITE ELEMENT METHODS FROM VARIATIONAL PRINCIPLES 
4.1 Introduction 

Much research has been done in the recent years to quantify 
uncertainties in engineering systems and their combined effect on the 
response. Theoretically, these uncertainties are modeled as random 
fields or random variables governed by joint probability density or 
distribution functions. In practice, the exact joint probability 
density functions are not always available; only the first few moments 
such as the mean, variance and correlations are known. The effect of 
these uncertainties on the system is ideally evaluated by examining the 
probabilistic character of the response, such as the probability of the 
response exceeding allowable limits. These limits are referred to as 
the failure surfaces. Recent research in reliability is focused on 
developing efficient techniques for this purpose. In general, these 
techniques involve much computation and are subject to various 
restrictions on the nature of the failure criteria. 

At the next level, estimates of response bounds, the level crossing 
rate or first passage time may be obtainable in some cases without 
extensive computations. At the easiest level of computation, the 
response statistics such as the mean, standard deviation and correlation 
coefficients are calculated [1,2]. These quantities are not only useful 
in themselves, but are also useful to calculate measures of reliability, 
e.g., the reliability index and reliability in terms of probability of 
survival [1-5]. In the past, the distribution functions for the 
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response were assumed and by using the response statistics, the 
probability of failure or survival was calculated [2]. Under a 
combination of uncertainties, the response distribution functions may be 
difficult to obtain. 

More recently, second moment reliability techniques have been 
formulated in the space of the input uncertainties. Here, the non- 
normal uncertainties are transformed to normal variables and the failure 
surface is described in terms of the response quantities. 

Transformation to normal variables is done to make use of the special 
characteristics of normal distributions, such as rotational symmetry and 
rapid exponential decay of the density function. If the failure 
function is linear in terms of the response, the probability of failure 
is expressed in a straightforward manner in terms of the response 
statistics [2,3]. If the failure function is nonlinear, approximating 
it by a quadratic function yields accurate reliability values [5]. 

For analysis purposes, time-invariant and time-variant 
uncertainties need to be distinguished. In the latter the probabilistic 
features, such as density functions and statistics, vary with time. For 
example, in a static analysis of structures, only time-invariant 
uncertainties, such as experimentally determined material properties, 
may be present* Conversely, in the dynamic analysis of structures, the 
forces such as earthquake excitations, wind and wave forces, and jet 
noise excitation of aircraft panels are often treated as time-variant 
uncertainties. Compared to time-invariant uncertainties, these are very 
difficult to quantify exactly, and assumptions such as Gaussian density 
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functions, stationarity, ergodicity and white noise characteristics are 
common [6,7]. In structural mechanics, the characterization of time- 
variant uncertainties in excitations, particularly seismic loadings, is 
an area of active research [8-10,14,15]. Expressed mathematically, 
these are differential equations with stochastic excitations. For 
linear systems the problem is tractable because of the applicability of 
spectral decomposition and superposition techniques [6,7,11]. For 
nonlinear systems, techniques such as equivalent linearization have been 
adopted with success [6,12,13]. A concurrent, albeit less extensive, 
area of research in structural mechanics is the one that deals with 
stochastic coefficients of both types. Problems with time-variant, 
stochastic coefficients have proven to be the toughest to analyze and 
still are an active research area. 

Apart from simulation techniques [2,16,17], a few non-statistical 
approaches are also currently available for solving problems with 
stochastic coefficients. These include numerical quadrature [18,19], 
second-moment analysis [20,21], the truncated hierarchy method [22], the 
method of moments [23], stochastic Green's function method [22], 
numerical solution of random integral equations [30,31] and the 
stochastic finite element method [29]. Merits and drawbacks of these 
and other methods are discussed in Refs. [24,26,29]. 

Recently, probabilistic finite element methods (PFEM) based on 
second-order perturbations have been formulated by the authors, for 
treating time-invariant, stochastic coefficients and excitations. These 
methods are based on second-moment techniques, so they are applicable 
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when the uncertainties are not too large and when the probabilistic 
density functions have decaying tails. Numerical results for various 
applications in trusses, bars, beams and plates have demonstrated the 
accuracy of PFEM, as compared with Monte Carlo simulation results [24- 
27]. These include both static and transient analyses of both linear 
and nonlinear structures, with random fields. Sample results for an 
elastoplastic cantilever beam with the uniaxial yield stress as a random 
field are given in Figs, la-ld [26]. The most appealing features about 
PFEM are its conceptual simplicity, ease of computer implementation, and 
computational efficiency. The special structure of the finite element 
equations, with features such as the symmetric stiffness matrices and 
the linear nature of the higher order equations, can be utilized to 
enhance the computational efficiency of the method for large-scale 
systems with many random variables. Compared to Monte Carlo 
simulations, computational requirements are often an order of magnitude 
smaller. 

This paper focuses on the development of efficient and accurate 
methods for calculating the response statistics in structural mechanics, 
making use of finite element modeling and solution techniques. By 
developing the PFEM equations from a variational principle, the 
randomness in the shape of the domain and boundary conditions can be 
treated. The PFEM equations are derived for linear continua from the 
potential energy variational principle in the next section. In Section 
4.3, the PFEM equations for nonlinear continua with large deformations 
are derived. It is shown that the final equations are similar to those 
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derived for the linear equations. Uncertainties arising from material 


and loads are accounted for in both formulations; randomness in the 
geometric properties such as shape is also included in the latter 
formulation. The PFEM equations are solved efficiently by numerical 
methods described in Section 4.4. The emphasis in this section is on 
numerical algorithms for computing the first and second order statistics 
of the displacement and stress and internal force. In Section 4.5, 
applications of PFEM to an elastoplastic plate with a hole and a turbine 
blade modeled by shell elements are studied. Results are summarized and 
discussed in section 4.6. 


78 


4.2 Development of PFEM for Linear Contloua from the Potential Energy 
Variational Principle 

The weak form which is obtained from the potential energy 
variational principle is: 


l '“(i.jAjktV.o* 1 'l iu i F i da - I su i h i dr ' 0 


(4.2.1) 


where the strain components are 


: ij U (i»j) 


A-f 1 3U < 3U 4 

d S f + — L'l 


in fl 


(4.2.2) 


and the stress components are given by the linear stress-strain law 


°ij * D ijk* e kZ 


in fl (4.2.3) 


The traction and prescribed displacement boundary conditions are given 
by 


o^n^ - h i on 33^ (4.2.4) 


and 


u 


i 


g i 


on 33 (4.2.5) 

g 


respectively. 5u^ is an arbitrary test function which satisfies 



6u i - 0 
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(4.2.6) 


on an 

S 

n is the domain, an^ and an^ are the traction and prescribed 
displacement surfaces, respectively, which satisfy 

an. n an - <p (4.2.7) 

h g 

In the above equations: u^ are the components of the displacement, Xj. 

are the spatial coordinates, nj are the components of the normal vector; 
^ijk& are the components of the material response tensor; F^, h^ and 
are the components of the body force, the prescribed traction and the 
prescribed displacement, respectively. Repeated indices denote sums and 
a comma denotes partial differentiation. 

The probabilistic potential energy variational principle (PPEVP) , 
which is a combination of the potential energy variational principle and 
the second-order perturbation method (l.e., the second moment analysis), 
embeds the probabilistic distributions, as reflected in the mean and 
covariance of the material properties, domain, boundary conditions and 
loading, to yield the corresponding means and covariances of the 
response in the variational statement. The basic idea of the second- 
order perturbation method in PPEVP is to expand each random function 
about the mean value of the random field b(x), denoted by T>(x) , and 
retain at most second order terms. That is, for a given small 
parameter 5, representing the scale of randomness in b(x) , the random 
function u^ is expanded about b via a second-order perturbation at a 
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given point x by 

<v 

o ' 2 ” 

u^ = + 5 ^ + 5 u £ (4.2.8) 

o ' 

where u^, u^ and u A are the zeroth, first and second order functions, 
respectively. Similar expansions are done for *i» h if n^ and g£. 

To simplify the subsequent development, the following abstract 
notations are introduced: 


<w,D,u> 

«w M 


[ W (i,j) D ijkf U (k,i) 


do 


(4.2.9) 


(w,P) - f w.P. do (4.2.10) 

~ ~ i i i 

a 

(w,h) - / w h dr (4.2.11) 

~ ~ 3 1 1 


Substituting the expanded functions into Eq. (4.2.1) and equating equal 
order terms, the zeroth, first and second order potential energy 
variational principles can be shown to be 


Zeroth Order 


o o 

<<5u,D ,u > 

/V Arf /V 


(5u,F°) + (5u,h°) a 

AJ /V /V# ^ £ 


(4.2.12) 


First Order (; terms) 



<Su,D°,u > « (5u,F ) + ($u,h ) - <5u,D ,u°> 

(V /v IV v /V v IV 3 v v v 
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(4.2.13) 


Second Order 



terms) 


<5u,D ,u > ■ (Su,F ) + (6u,h ) - <6u,D ,u > - <6u,D ,u > (4.2.14) 

v v v v v v v v v v v v v 


Remark 2.1 The arbitrary test function 5u satisfies 5u * 0 on the 

— — V V V 

boundary 3R . 

8 

Remark 2.2 All the functions with a superscript " o " are deterministic 
functions (l.e., evaluated at b) whereas functions with superscripts 
and are random functions characterized by the random field b(x). 


Remark 2.3 Equation (4.2.12) is the standard deterministic variational 
statement and therefore, the usual Galerkln finite element procedure can 
be employed directly. Once u° is determined from Eq. (4.2.12), the 

I " 

random functions u and u can be determined using Eqs. (4.2.13) and 

V /V 

(4.2.14), sequentially. 

»i» » 

It should be noted that the random functions 0 , F , h and £ and 
the functions with the superscript are, in general, described 
through spatial expectation and autocovariance functions. Therefore, In 
addition to the usual finite element approximation of the displacement 
field, the random field is also discretized with q shape functions. To 
Incorporate the spatial expectation and autocovariance functions into 
the formulation, the discretized random variables b K are expanded 
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about t> via second moment methods. The first order variation at each 
& 

point x is denoted by db„ - £Ab - g(b -b ). 

~K ^ K K. K 

For consistency with the finite element approximation and to assure 
the accuracy of the approximation of the random field, the random 
functions D, F, h and g; which are, in general, functions of b(x) and 

m /w /w /v 

x, are first discretized with the same q shape functions. For example, 
the finite element approximation of D, (the coefficient 1/2 has been 
added to the quadratic term (Eq. (4.2.18)) so that it is consistent with 
conventional descriptions of second moment analysis) is given by: 

D « D° + + £ 2 D (4.2.15) 

M O# <-<W 


or 


D * E # t (x){D° + sd’ + 5 2 DJ (4.2.16) 

** I ~ ~ I ~I J 

O th t 

where D T denotes the I nodal value of D evaluated at b, D_ denotes the 

fi 

first order variation of D(x_,b) due to variations Ab„, D_ denotes the 

second order variation and $ (x) are the q shape functions. The random 

l ^ 

i 

variables are then expanded in terms of the random variables b^ by 


f 



q 

z 

K-l 


( 2i>k Ab K 


Si 


q 

z 

K,L» 


( 2i^KL Ab K Ab L 


(4.2.17) 


(4.2.18) 
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f " 

respectively. The nodal values of (D_)„ and (D_)™ can be obtained by 

K ~I KL 

partial differentiation or by a least-square fit. Similar definitions 
hold for P, h and g. 

The displacement is discretized similarly, however, with NUMEL 
elements and NUMNP nodes with each node having NDOF degrees of freedom, 
via (c.f. Eq. (4.2.8)), where 


u 


NUMNP 

Z 

A-l 

NUMNP 

Z 

A-l 

NUMNP 

Z 

A-l 


N.(x)d° 
A ~ ~A 

N,(x)d" 
A ~ ~A 


(4.2.19) 

(4.2.20) 

(4.2.21) 


and N^ are the displacement shape functions. The first and second order 

I •• 

variations of d A and d k are defined by 
~A ~A 





q 

z 

K-l 


( ~A } K Ab K 


(4.2.22) 


and 


~A 


q 

Z 

K,L- 


( VKL A V b L 


(4.2.23) 


respectively. 

Substituting the above Galerkln/f inite element approximations into 
the zeroth, first and second order variational statements, the finite 
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element matrix equations can be obtained using Eqs. (4.2.19-21) and the 
arbitrariness of 6d, 

6d m (6d^,6d2» • « • (4.2.24) 

It should be noted that only NEQ test functions are required since no 
test function is needed on the parametrized boundary 3R ; NEQ is the 
number of finite element equations to be solved. 

Zeroth Order 




(4.2.25) 


where 

K - {K^} - <N a , D°, N b > (4.2.26) 

f° - { f°> - (n a , F°) + (n a . h°) 3 

- <N a , D°, (4.2.27) 

d° - {d°} (4.2.28) 


where the subscript A takes the values of 1 to NEQ, B is summed from 1 
to NEQ and C is summed from NEQ + 1 to NED; NED is equal to NDOF times 
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NUMNP. 


First Order (for each Ab^) K ■ 1, . .., q 


r 


K d 

/V /V 


K 



(4.2.29) 


where 


i 



< f A>K 



~K J + (N A* ~K J 3 " <N A’ 



- <» A . B 9 * V ( 8c>k 


(4.2.30) 


and 


4c ' < d B>K <4 - 2 ' 31) 


Second Order (K and L are summed from 1 to q) 


K d - f 


(4.2.32) 


where 


{f A } KL Ab K Ab L 


1 2 (N A> ~KL^ + 2 (N A» ~KL^3 “ 2 <N A’ Hrl* ~ > 



< N A* 5k» “ 2 <N A’ s » N C > ^ g C^KL^ Ab K Ab L 
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(4.2.33) 


and 


i ■ 1 < Via 4 V b L 


(4.2.34) 


The mean and autocovariance matrices for the nodal displacement are 
defined by 


E[d] - / d(b)f( 5 ,b)db 

^ * /V M A# 


(4.2.35) 


and 


C0V(d tA’ d jB ) (d U ~ d lA )(d jB - d j»> £( 5- b > db (4 - 2 - 36 > 

respectively, f is the joint probability density function which is 
dependent on 5 . b is the random-variable vector 


b - (bj, b 2 , ..., b ) 


(4.2.37) 


With Eq. (4.2.35), the second order accurate mean value of d is shown to 
be 


E[d] « d° + d 

rs* 


(4.2.38) 
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where d Is given by 


d 


1 

2 


q 

Z 

K,L-1 


S*1 Co '' (b K' b L > 


(4.2.39) 


From Eq. (4.2.36) the first-order autocovariance of the displacement can 
be shown to be 



To simplify the computational procedure of Eqs. (4.2.32-34), 
integrate Eq. (4.2.32) over the range of the random variables so that it 
is replaced by 



(4.2-.41) 


where d is given in Eqs. (4.2.39) and f is given by the right-hand 
side of Eq. (4.2.33) with Ab^Ab^ replaced by Cov(b^,b| j ). 
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4.3 PFEM for Nonlinear Continua with Large Deformations 

In this section, the PFEM formulation for large deformations of 
hyperelastic materials in bodies of random shape is considered. The 
simplest form of the weak form for nonlinear elasticity with large 
deformations is s 

/ 5e T a dn - / 5u T F dft + / 5u T h dn (4.3.1) 

ft n an h 

where e, a, F and h are the nonsymmetric measures of the strain, first 

/v /w 

Piola Kirchhoff stress, body force and traction, respectively; ft and 

3ft. are the domain and natural boundary, respectively, in the initial 
n 

configuration. Furthermore, the strain measure is given by 

e - Vu - G - I (4.3.2) 

where G is the deformation gradient and I is the identity matrix. The 
stress for a hyperelastic material is given by 

a - 4(G) (4.3.3) 

/V /V (V 

where 


t 


3W 

3G 


(4.3.4) 


and W is the strain energy density function. Randomness in material. 
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geometric properties and loads are represented by the random field 
b(x). To Incorporate the random domain and boundaries Into the 
formulation, Eq. (4.3.1) can be rewritten as 

f de^c J dR ■ / du^F J dR + / du^h J dA (4.3.5a) 

J /v «v y J » s# y J m g 

R R 3R. 

employing the following mappings of the original domain 3 and boundary 

33, onto the reference domain R and boundary A: 
a 


da - J dR 
v 


(4.3.5b) 


and 


dr - J dA (4.3.5c) 

8 

The displacement is approximated by a second order perturbation about 
the mean random field b at a given x: 

o f 9 " 

u a u + £u + £ u (4.3.6a) 


where f represents the scale of randomness in b(x). Similarly, e, F, 
h, and J g are expressed as second order perturbations. The stress is 
expanded as 


o ' o » 2 " » • o " 

o * * + ?(<(» + Ce)+§(]1» + Ce +Ce) 

IV /w ’ /V /w /v /v A# 


(4.3.6b) 
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where the first elasticity tensor is given by 


C 


3 G3 G 


(4;3.6c) 


The derivatives of nonlinear functions such as 4) are given in Ref. (28 J; 
however for nonlinear elasticity, the derivatives can be worked out 
explicitly. In addition, since the boundaries are parametrized, it is 
necessary to perturb the virtual strain in Eq. (4.3.5a) as this is a 
function of the domain geometry. Thus, 


o * ? " 

Se = Se + g6e + g 5e (4.3.6d) 

The randomness in domain and boundary geometry is taken care of by the 
Jacoblans J y and J g , respectively, in Eq. (4.3.5a). 

Substituting Eqs. (4.3.6a) through (4.3.6d) in Eq. (4.3.5a), the 
PFEM equations are obtained: 


Zeroth Order 


/ 5e°^i|» 0 J 0 dR * / 6u^F°J° dR + / 6u^h°J° dA (4.3.7) 

* A# /W V 7 A< A# V J M A« S 

R R 3R h 

First Order 


R 



Second Order 


/[*£ V J v + 6 £ 1( 4 + S°£ )J v + 5 £ Vi 


o \,o . . *T.o ' 


R 


+ 6e° T (t|> + cV + C° e )J° + 6e° T (^ ' + C 0 e')j' + 6e° T d»°vJ IdR 

iw m m ^ Y <v /w V ^ A* V J 


/[«sY< 

R 


X f f T O ,# n 

Ju F J + fiu F J IdR 

v m v ^ 


+ / [Su A h J° + 6u A h J + 6u^h°J IdA 

4 l /V m S ** ts* 3 'v <v O J 


(4.3.9) 


It is noted that the test function 5u satisfies 6u » 0 on the 

<v <v 

parametrized prescribed displacement boundary 3R and the intersection 

& 

of the boundaries 3R^ and 3R^ is a null set* 

Since Eq. (4.3.7) is the deterministic virtual work principle, 
standard techniques such as Newton Raphson iteration can be employed for 
its solution. After determining u°, e° and £° from Eq. (4.3.7), the 

f I f »* « M 

random functions u , e , ib and u , e , i|r can be determined from Eqs. 
(4.3.8) and (4.3.9), sequentially. It should be noted that the random 

• •tit i 

functions * , C , F , h , J and J and the functions with the 
x ’ 1 « ’ I. ’ s v 

superscript are, in general, given and described through the spatial 
expectation and autocovariance functions. Similar to the. previous 
section, the random functions are discretized with q shape functions. 

For example, the finite element approximation of the first 



elasticity tensor, C, is given by: 
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(4.3.10a) 


or 


C * 


E # I (x){C I + 
1-1 



(4.3.10b) 


where C° denotes the I c ^ 
~I 


nodal value of C evaluated at b. 


si - j. <Si>M 4b M 


(4.3.10c) 


and 


C - -r 


q 

z 

M,N- 


( £lW Ab M Ab N 


(4.3.10d) 


respectively. Similar definitions hold for F, h, jji, J y and J g . 

Similarly, the displacement field is discretized, however with 
NUMEL elements and NUMNP nodes with each node having NDOF degrees of 
freedom (see Eqs. (4.2.8) and (4.2.19-23)). 

The elemental strain is expressed as 


e - B d (4.3.11) 


where B is the discretized gradient operator. The strain perturbation 
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is: 


o 

e « e 



2 " 
+ 5 e 


(4.3.12) 


or 


e * B°d° + E(B'd° + B°d') + £ 2 (B d° + bV + B°d ) (4.3.13a) 

fS* /NS /V * /NS /NS M SW ^ /V /NS (NS IV AS 

and from Eq. (4.4.3.6d) 

5e » B°6d° + eb' 6d° + E 2 B 5d° (4.3.13b) 

IV /V IV /NS /NS /V /NS 

Substituting Eqs. (4.3.10) through (4.3.13) in Eqs. (4.3.7), (4.3.8) and 
(4.3.9) and using the arbitrariness of 6d° the zeroth, first and second 

/NS 

order equilibrium conditions are obtained, respectively. 

Zeroth Order 

/ B° T 4°J 0 dR - / N T F°J° dR + f N T h°J° dA (4.3.14) 

* /V A* V / /V AS V S /NS /NS Q 

R R 3R h 

First Order (for each b^, M * 1, ..., q) 

K <L, - f’ (4.3.15a) 

INS /Nsfl /Nsfl 


where 



94 


f 



/ R /[£>° + £ 0 <V M ] dRt 


1 /[yj * 

3R h 


r _’T o _0 e „°T, ’ 

- / 5 4 J v dR " / 5 

R R 


+ C 0 B’d°)J 0 dR - f B° T U; 0 ( J * ) u dR(4.3.15b) 

~ ~ ~ V ~ ** V M 


with 


C° - D° + T° (4.3.15c) 

*** *+ ** 

resulting In 

K - K_ + K_ (4.3.15d) 


where 


„ , _oT„o_o T o 

1L ■ ] B DBJ dR 

/M U J A* Cv V 

R 


and 


(4.3. 15e) 


K - / B 0T T°B 0 J° dR (4.3. 15f) 

~G J ~ ~ ~ v 

R 

where D° and T° represent the material response and Initial stress 

A* 

matrices respectively. Explicit expressions of these matrices for a QB1 
element are given in Ref. [38]. 



Second Order (M and N are summed from 1 to q) 
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K 

where 


and 


Eqs. (4 



<4. 


£ “ {/ H C 2 ~MN J v + + 2 ~ ^ J v^MN^ dR 

K 


*1 S T (i !Ws + VVn + 1 5°< V*]<* 

3E h 

- / .; t % * s°*d° + sYi>: <** - / &y<v s “ 

K. R 


„oT, 1 


• » « • « » 


1 S ”‘<2 1 «n + £m5Y * * 2 £°W° + 

R 


/2 01 <*m + £V + £°A»Vn dR 

R 


- / 7 »° T t^Vai dR - / 7 (4 ‘ 

R R 


~ “ 2 2MN^ b M^ b N 


(4, 


.3.14) through (4.3.16) are solved in sequence. 


3.16a) 


dR 


3.16b) 


3.16c) 


The second order accurate mean value of d is 
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E[d] = d° + d 

A/ (V 

(4.3.17a) 


«• 

where T is given by (c.f. the previous section) 



~ " *2 M N-l C ° V(bM * l>N) 

(4.3.17b) 


and the £irst order accurate autocovariance of the displacement is 


C ° V ( d iA ,d jB^ “ * ^ d iA^M^ d jB^N C ° V ^ b M* b N^ 

(4.3.17c) 


The mean and autocovariance matrices for the element 

strains can be 


obtained from the following: 



H 

E[e] ■ e° + e 

M# M M 

(4.3.18a) 


where 



e * B d 

A# /V 

(4.3.18b) 


and s is defined by 

/■w 



~ * M n 1 ^ + + 1 2 d MtJ Cov(b M ,b N^ 

M, N»1 

(4.3.18c) 


Thus, 



Cov <Si»£j> * ^ ^~i^~j^N ( ' ov ^ b M ,b N^ 

M,N»1 

(4.3.18d) 
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with the definition 


' .o .o, ' 


^~i^M " ^i^i + ~i^~i 


i * 1, NUMEL, no sum on i (4.3. 18e) 

M a 1| • i • | 


where subscript "i" denotes the i tb element and NUMEL is the total 
number of elements. Similarly, the mean first Piola-Kirchof f stress is: 


E[o] ■ o° + o 


(4.3.19a) 


where 


o o 

2 m t 


(4.3.19b) 


and 


Z - l ,1 ti hrn + + * i STwr 

M.N*! 


1 „o o 


+ S + 2 £ 5 it*] Cov(b H’ b H ) 


(4.3.19c) 


and the stress autocovariance is defined by 


Cov ^2i»2j^ * 2 ^2i^M^Sj^N Cov ^ b M* b N^ 

M,N*1 


(4.3.19d) 


with the definition 



f 


( 2i>M 


+ 
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i ■ 1 , NUMEL , no sum on i 

M ■ 1 f • * • j (j 


(4.3.19e) 
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4.4 Computational Aspects of PFEM 

4.4.1 Transformation of the Covariance Matrix 

In linear problems, it can be observed from Bqs. (4.2.25), (4.2.29) 

o * — " 

and (4.2.41) that the determination of d , d,, and d involve one 

~ ~K ~ 

factorization of K and q + 2 forward-reductions and back-substitutions. 

The latter operations consist of evaluation of d° in Gq. (4.2.25), q 
» 

evaluations of d v in Eq. (4.2.29) and one more evaluation of d in Eq. 
(4.2.41). Even though the above computations are compatible with the 
elemental discretization and nodal assembly procedures that characterize 
finite element theory and software, the number of matrix multiplications 
is proportional to q(q + l)/2. This would be unacceptably expensive. 

The large number of computations arise from the double summations in 1 
and j in Gq. (4.2.33). To remedy this situation, the covariance matrix 
Cov(b^,bj) is transformed to a diagonal variance matrix Var(c i ,Cj) such 
that 


Var(c i ,Cj) - 0 


for i * j 


(4.4. la) 


and 


Var(c i ,Cj) - Var(c i ) 


for i » j 


(4.4.1b) 


Therefore, the number of matrix multiplications is now proportional to 
q. The above involves the solution of the eigenproblem 



G t(> ■ i)>A 
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(4. 4. 1.2) 


where G^j and denote Cov(b^,bj) and Var(c^,cj), respectively, 

and ip is a constant transformation matrix with the following properties: 

j m 1 (4.4.1.3a) 

A - (4.4.1.3b) 


and 


b ■ Ac or c 

/Nrf (WV A# 



(4. 4. 1.4) 


^ is the q x q identity matrix and c is the transformed q x 1 random 
variables vector. The mean and variance are 


Etc] - d» T E [b J 

IV IV IV 


(4.4. 1.5) 


and 


Var(c) =* diagonal terms of A (4.4. 1.6) 

*** IV 

With Eqs. (4.4.1.3b) and (4. 4. 1.4), the mixed derivatives appearing 
in Eq. (4.2.40) reduce to second derivatives, and the covariance matrix 
is replaced by the diagonal variance matrix in Eqs. (4.2.33) and 



101 


(4.2.40). 

Analogous to modal analysis In structural dynamics problems, only a 
few modes n(n < q) (i.e., Var(c i )) are required to capture the major 
characteristics of the probabilistic distribution. However, the highest 
eigenvalues have to be employed. This is in contrast to the modal 
structural dynamics problem, wherein the lowest eigenvalues are used. 

As PFEM Involves, essentially, a set of sensitivity equations with 
respect to c, recent techniques in design sensitivity analysis can be 
adapted easily. One such technique is the adjoint method in mechanical 
design [32-35]. In this method, the first and second order derivatives 
of the objective functions and constraints are calculated w.r.t. the 
design parameters, with minimal computations of the first and second 
order equations. 

4.4.2 Adjoint Method in PFEM 

Consider a typical function dj(c,d) involving the displacements d 

T (V 

and the random variables c. Chain differentiation yields 

/N* 

T 1 

[d) ] * iji i * 1, ...,n (4.4.2.1a) 

yj c J y c. y d ~c. 9 9 

1 l ~ l 

where the subscript denotes the derivative with respect to c^, and 
T 

iK ■ (i|>. > ..., ip, , ..., iJ/j ) (4.4.2.1b) 

5 a l a k NEQ 

Substituting Eq. (4.2.29) in Eq. (4.4.2.1b), the explicit equation 
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[4» K _1 f i - 1, ..., n (4. 4. 2. 2) 

c. c. a ~ ~c. 

i i ~ i 

is obtained# Usually, in the direct method, the above equation is 
evaluated for each random variable c^, involving f n f solutions of the 
linear equation (4#4#2#2). In the adjoint method, X is selected to 
satisfy 


K \ - ip, (4. 4. 2. 3) 

/v Q 

Then, Eq. (4. 4. 2. 2) can be rewritten as 

[*l c * ♦ «. + ^ £e 1 ■ 1» •••» a (4. 4. 2. 4) 

C i C i C 1 

The adjoint problem, Eq. (4*4. 2. 3), is solved only once in this 

method. In the direct method, 'n' solutions of Eq. (4.2.29) are 

required. This is the advantage of the adjoint method over the direct 

» 

method# Both methods require 'n 1 inner products with f , in Eqs. 

~°i 

(4.4.2. 1) and (4. 4. 2. 4), respectively. However, it has been shown that 
when the number of functions iji is more than the number of random 
variables, the computational advantage of the adjoint method is lost 
[33,34]. By solving 'n' adjoint problems, the second order sensitivites 
can also be evaluated [33,34]. It should be noted that the adjoint 
method is applicable to nonlinear problems as well, as the first and 
second order equations are still linear. 
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4.4.3 Displacement Derivatives in PFEM 

In mechanics, one is often interested in the response in only a 

portion of the entire domain. Stress concentrations, plastic flow and 

strain localization effects are some examples. Similarly, in 

probabilistic analysis, one is interested in the probability of failure, 

which usually initiates in a small domain. This translates into a few 

nodal displacements and element strains and stresses. In such cases, 

the adjoint methodology can be used to reduce the computations in PFEM 

equations for the evaluation of derivatives. 

The adjoint method can be used to calculate the displacement 

derivatives of the k 6 * 1 component of the displacement vector d, denoted 

** 

(k) 

by d . This is done by substituting 

tji - d (k) (4.4.3. 1) 


in Eq. (4. 4. 2. 4). Thus, 


M< k >] - d <k > ♦ x 1 f' 

C i C i ~ ~ C i 


(4.4. 3. 2) 


where X is obtained from the adjoint problem 


K X - d 


(k) 


(4. 4. 3. 3) 


Interestingly, the right hand side of Eq. (4.4. 3. 3) is a Boolean vector, 
with unit value at the k c ^ component. Therefore, the adjoint problem 
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for displacement derivatives can be Interpreted as a linear structural 
problem with the same tangent stiffness and a unit load at the k ttl 
position of the external force vector. The displacement covariance Is 
then obtained from Eq. (4.2.39), where the derivatives are with respect 
to c and the covariance matrix Is replaced by the diagonal variance 
matrix (Eq. (4. 4. 1.6)). 

•t 

In the direct method, the second-order term d in Eq. (4.2.41) is 

Al 

obtained by a single solution, irrespective of the number of random 
variables c^. This Is because of the summation of the second-order 
displacement derivatives In Eq. (4.2.39). 

In comparison, the adjoint method may require more computations to 

M 

compute the second-order term d . This term requires the first 

/V 

derivatives of the displacements, over the entire domain . In such a 
scenario, the adjoint problem (cf. Eq. (4. 4. 3. 3) has to be solved for 
each component of the displacement vector d, resulting in more 

/v 

computations. If the size of the vector d is small when compared with 
the number of random variables c^ f the adjoint method will require fewer 
computations than the direct method. Thus, the selection of the adjoint 
method over the direct method depends on (i) the number of displacement 
components considered, (ii) the number of random variables c^, i * 1, 
..., n and (iii) the size of the displacement vector d. 

/v 

It Is to be noted that the adjoint problem is always linear, 
irrespective of the primary problem. It has been noticed that the 
second-order term contributes very little to the mean displacement 
calculations [24-26]. If the second-order term is neglected, then the 
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adjoint method for the first-order mean and covariance would involve 

solutions of only two equations viz., Eqs. (4.2.29) and (4.4.3. 3) and 

? 

'n' inner products with f in Eq. (4. 4. 3. 2). The adjoint PFEM method 

~ C i 

for displacements is applicable to linear and nonlinear materials, with 
the use of stiffness and tangent stiffness matrices, respectively* The 
first order mean and covariance of displacements d^) f k, * 1, • •*, N 
where N is the number of displacements of interest, are 

E[d (k) ] - d° (k) k - 1, ...» N (4.4.3.4a) 

Cov(d (k) ,d (A) ) - { Z [d (k) ] [d (i) ] Var(c )} (4.4.3.4b) 

1 . c c r J 

r-1 r r 

4.4.4 Stress Derivatives in PFEM 

The first derivative of the stresses with respect to the 
probabilistic variables, in any element, can be expressed in terms of 
the displacement derivatives in this element which are first calculated 
by the adjoint method. For a four node, 2D continuum element this 
requires the solution of eight adjoint problems. For linear materials, 
the stress derivatives with respect to the transformed vector c, in a 
given element are: 


2 Cj 


d«d 


„ _0 ,0 , o o 

C B d + C B d 

/N* M #v(J 


(4.4.4. la) 


resulting in 



[?] 


o o 

^ + CBd 
o ~ ~ ~c , 
d=d i 


106 

(4.4.4.1b) 


In the absence of random material properties, the first term in Eq. 
(4.4.4.1b) drops out and in the absence of random geometric properties, 
the second term drops out. In the case where only loads are random, 
both terms drop out. 

For nonlinear materials, Eq. (4.4.4.1a) cannot be easily 
evaluated. One strategy is to replace these derivatives by their 
finite-difference counterparts [24] 


o' j * - (o° + - o°") I 

~ C i|d-d° ^i ~ ~ |d-d° 


(4.4.4.2a) 


with the definitions 


0+ / ° , . X 

a * o(c + Ac. ) 

<"W M 


(4.4.4.2b) 


0° ■ o(c° - Ac ) 

a* /v <w 


(4.4.4.2c) 


and 


A £i 


(o, 


Ac 


i* 


0)' 


(4. 4. 4. 2d) 


The derivatives of the tangent constitutive matrix in Eqs. (4.3.16b) can 
also be approximated similarly. For a general nonlinear material, Eqs. 
(4.4.4.2b) and (4.4.4.2c) have to be evaluated by solution of the zeroth 
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order equation, with the appropriate values of c* However, for 
elastoplastlc materials with random material properties, these 
quantities can be evaluated in the course of the zeroth order solution 
by the radial return algorithm [36]. Additional arrays to store the 
stresses in Eqs. (4. 4. 4. 2) are all that is needed to achieve this. 
Essentially, the functional relationships 


o+ 

~t+At 


- o+ .o . 
2<£t ’ St’ A £t ) 


(4.4.4.3a) 


and 


~t+At 



(4.4.4.3b) 


hold, for each c^ [37]. The subscripts 't' and ’t+At' refer to two 
successive time steps, in the evolution of the stress history. The 
first order mean stress and covariance of stress are expressed as 


E [a] 


o 

“ a 


(4.4*4. 4a) 


and 


Gov(a i ,o^) - { Z [o i ] r [cr* ] r Var(c f )} (4.4.4.4b) 

r*l 

where <j* end are any two components of the elemental stress 


vector <j 
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4.5 Applications 

The usefulness of PFEM is demonstrated here by three 
applications. In the first two cases, the problem studied is an 
elastic-plastic plane continuum with a circular hole. In the first 
application (Fig. 2), the uniaxial yield stress and the uniform 
compressive load are assumed to be two independent stationary random 
fields, with an exponentially decaying correlation function. The load 
is discretized into 12 components, applied at the nodes, both for 
deterministic and random analyses. The yield stress is assumed to be 
radially correlated in an exponential manner (Fig. 2). The domain of 
the plate is divided into 15 ring-like bands and in these bands the mean 
yield stress is assumed constant in time and space. This results in 15 
discretized random variables for the yield stress. The load is quasi- 
static and linearly increasing with time, with 8 load steps. The 
displacement and stress statistics are studied in each of these load 
steps. 

The mean and variance of the compressive stress at load step 8, 
along the x-axis, are plotted in Figs. 3a and 3b, respectively. The 
elements near the hole are plastic at this load, as can be seen from 
Fig. 3a. The variance of the stress is maximum, therefore, near the 
hole and this is seen in Fig. 3b. These results compare very well with 
the Monte Carlo simulation (MCS) results. The coefficient of variance 
is 10% for the stress. The mean and variance of the stress near the 
hole (Point B) are plotted in Figs. 4a and 4b, respectively. The mean 
stress is plastic after the second load step. Thereafter, the mean 
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stress Is almost constant and is equal to the uniaxial yield stress 
value, except for a slight hardening effect. The variance of stress 
rises rapidly when the yielding occurs (i.e., elastoplastic state of 
stress) and thereafter it rises gradually. The maximum coefficient of 
variation for the stress at Point B is 10Z. The stress correlation 
along the x-axls, w.r.t. the stress at Point C, is plotted in Fig. 4c. 
Interestingly, the correlation is almost zero near the hole and there is 
inverse correlation near the fixed end of the plate (cf. Fig. 2). This 
suggests that the stresses in the plastic state are not correlated with 
the stress near Point C. This can be explained by the fact that beyond 
yielding the stress remains practically constant for this material 
because E/E^ * 100. Near the fixed end, the stress is almost zero and 
the negative correlation implies that this stress will statistically 
increase when the stress near point C decreases. By studying the 
effects of the two random fields separately, it is noticed that the 
random load effect in terms of the variance, is spread wider over the 
elements along the x-axis than the random material effect, and the 
effect is mainly near the hole. As is to be expected, the variance of 
the stress under the combined effect of the random material and random 
load is additive. Similarly, for the mean values of stress, the second 
order effect is additive. 

The second example concerns the cyclic loading of the same plate 
with only the yield stress as the random field. Mechanical and 
aerospace components are usually subjected to thousands of cycles of 
stress, resulting in fatigue. The material properties usually show some 
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degradation with time In these components. A modest attempt Is made 
here to see If there is a large variation of the response statistics 
after 3 cycles of loading and unloading. The mean and variance of the 
displacement at node A is plotted as functions of the load step in Fig. 
3a. The mean displacement is sinusoidal, resembling closely the forcing 
function. The variance of the displacement is zero until the plate 
begins to yield in compression. After this, the variance jumps to a 
higher value and remains steady until the yielding in tension begins. 

At this point, there is a sharp drop in the displacement variance and 
after that the variance stays at a constant level. This phenomenon 
repeats every cycle. There seems to be a gradual buildup of the value 
of the displacement variance during every cycle, particularly under 
compression. The maximum coefficient of variation in these cycles is 
2Z. 

The mean and variance of the stress at Point B are plotted in Fig. 
5b. The mean stress is periodic, with a slight flattening at the top 
and bottom. This flat region corresponds to the plate yielding, in the 
mean sense. The variance of stress at Point B is periodic, and behaves 
similar to the displacement variance. The coefficient of variation is 
10Z. 

The stress variance exhibits spikes whenever yielding is about to 
commence. The variance drops to a near zero level in these downward 
spikes. This phenomenon can be explained from the elastic-plastic 
behavior of the plate under stress reversal. To do this, three 
deterministic solutions of the stress in the plate are shown in Figs. 6a 
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and 6b at various loadsteps, under cyclic loading. These correspond to 
yield stress (o y ) values of 26,000, 25,000 and 24,000 psi, respectively. 
During the first yielding, in compression, the magnitude of the stress 
for Oy « 26,000 is maximum and this stress plot, in Fig. 6a, lies 
outermost. Before the next yielding in tension commences, there is a 
crossover of the three curves. This crossover is necessary because the 
magnitude of the stress at yield, for the highest yield stress value, is 
always the highest. The crossover repeats twice for each cycle of 
loading. This translates into a very small variance of stress near the 
crossover regions because the variation of stress w.r.t. the yield 
stress is near zero. The spikes in the displacement variance can also 
be explained similarly. 

The third application studied is a turbine blade with random load 
along the edge, random yield stress and random length of the blade. The 
problem statement, along with the details of the random fields, are 
given in Fig. 7. The load is quasi-static and linearly increasing, with 
15 load steps. The expectation and deviation of the displacement are 
plotted in Fig. 8a; the coefficient of variation is plotted in Fig. 

8b. It is noticed that the first two steps are elastic and beyond that 
the blade starts yielding. Due to this yielding, the expectation is 
nearly flat beyond the second load step. In the elastic region, the 
maximum contribution for the deviation comes from the random load 
followed by the random length. In the elastoplastic region and beyond, 
the random yield stress affects the deviation most. The combined 
deviation has a maximum coefficient of variation of 13% and this occurs 
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just after the initiation of the yielding. 

The stress statistics are plotted in Figs. 9a and 9b. The stress 
deviation is largely due to the random load in the elastic region, as in 
the case of displacement deviation; the effect of the random length is 
very small in both the elastic and plastic regions; the random yield 
stress causes the most deviation in the plastic region. The maximum 
coefficient of variation is 8%, at the last load step. 
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4.6 Conclusions 

The PFEM techniques for linear and nonlinear materials, Including 
elastoplastic materials, yield efficient and reliable statistics of the 
response quantities of interest. The direct solution of PFEM equations 
may require a substantial number of computations for large systems. By 
making use of features such as the eigenvalue orthogonalizatlon and 
selection of only a few highest eigenvalues, the adjoint methodology, 
and superposition of random fields, these computations can be 
drastically reduced. The results obtained here show that the first and 
second order variances in response which are obtained by this form of 
the moment method agree well with Monte Carlo simulations when material 
properties such as the yield stress are random variables. The results 
also seem to suggest an increase in the response variance even with a 
small but steady degradation of material properties after several 
cycles. However, this needs to be investigated further. As discussed 
in the introduction, based on the response statistics, reliability 
measures can be calculated. In addition, the response gradients with 
respect to the random variables are calculated in the course of PFEM 
calculations. These are also useful in reliability calculations and 
probabilistic design. 
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CHAPTER 5 


APPLICATIONS OF PROBABILISTIC FINITE ELEMENT METHODS 
IN ELASTIC/PLASTIC DYNAMICS 

5.1 Introduction 

Design methods for engineering problems are, in general, based on 
deterministic parameters. In practice there are often uncertainties 
associated with parameters such as: material and geometric properties, 

forces, and boundary conditions. Although, in most situations the 
uncertainties may be small, the combination of these can lead to large 
and unexpected excursions of the response, particularly in multi- 
component systems. In the context of failure and reliability analysis, 
this phenomenon is of obvious significance. In the past, problems with 
uncertainties have been studied to provide an Insight of the statistical 
response variations, with methods like sampling [1-4], numerical 
integration [5,6], second-moment analysis [6,8] and stochastic finite 
element methods [6,9-12]. The choice of the appropriate method depends 
on the nature of the problem and this was briefly discussed by the 
authors in Refs. [12,13]. Typically, the uncertainties are modelled as 
random quantities governed by probability density functions, and that is 
also the case here. 

A survey of the existing literature shows that, with the exception 
of the methods based on sampling, the other methods are limited to 
linear problems. Moreoever, techniques for handling random fields, 
where the randomness is spaced over the contlnua, are even scarcer. The 
authors have recently extended probabilistic finite element methods PFEM 

[13] to linear and nonlinear continua in both static and transient 
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settings. 

A schematic of the PFEM is presented in Fig. la. In the PFEM [13], 
the random fields, characterized by the mean, variance and 
autocorrelation functions, are discretized to obtain the mean vector and 
the covariance matrix. For a correlated random field, the covariance 
matrix will be a full matrix and therefore it will require too many 
computations. To remedy this, the correlated vector is transformed to 
an uncorrelated vector by an eigenvalue orthogonalization procedure 
resulting in a diagonal covariance matrix, and therefore, fewer 
computations. This transformation procedure gives rise to a set of 
modes and corresponding eigenvalues. It is shown that only a few of 
these modes are sufficient to obtain a converged PFEM solution. 

Finally, the PFEM involves solution of a set of deterministic FEM 
equations to obtain the mean, variance and autocorrelation of the 
response. 

In this paper, two applications of the PFEM in elastic/plastic 
dynamics with random material properties are studied in detail. The 
discretization of the random field depends on factors such as the 
inhomogeneity of the randomness and the extent of the spatial 
correlation. The necessary guidelines for the discretization are 
discussed in the next section. In Section 3.3, the choice of the number 
of modes necessary for a converged PFEM solution is discussed. The 
computational efficiency and accuracy of this method are compared with 
those of Monte Carlo simulation with a first-order filter [6,7] in 
Section 5.4, along with the conclusions. 
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5.2 Random Field Discretization 

Let b(x) represent the random field. In PFEM, b(x) is approximated 
by discretization as 


q 

b(x) - I N. (x)b. (5.2.1) 

~ 1-1 i ~ 1 


where N. (x) represent the shape functions and b, the discretized values 
I ~ 1 

of b(x) at x^, i - 1, ..., q. It follows from Eq. (5.2.1) that 


q 

db(x) - E N (x)db (5.2.2) 

~ i-1 1 ~ 1 

2 q 

db (x) - E N. (x)N. (x)db. db (5.2.3) 

i,j-l 1 ~ j ~ 1 J 


where 


db - b - b 
i i i 


(5.2.4) 


and b^ represent the mean values of b^ (also denoted by the expectation 
operator £[•])• From Eq. (2.1) the expectation and the covariance of 
b(x) are, by definition. 


E[b(x)] - / b(x)f(b)db 

/W * M MW 


(5.2.5) 


E N.(x)E[b ] 
i-1 1 ~ 1 


(5.2.6) 
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and 


Cov(b(x^) ,b(x^)) - / (b(x fc ) - b(X| t ))(b(x^)-b(x^) )f (b)db (5.2.7) 


* » J“1 


(5.2.8) 


where f(b) is the multivariate probability density function; x. and x 

~ ~K 

are any two points in the domain of x. 

From second-moment analysis [6,8], the mean of any function 
S(b(x) ,x) at any point x, , and the covariance of the function between 

/v /v /M 

any two points x. and x can be written as 


E[ Sk ] 



q 

E 


i.j-1 






Cov^.bj) 


(5.2.9) 


and 


Cov(S k ,S £ ) 



(5.2.10) 


where 


S k - S(b(x),x k ) (5.2.11) 

and the superposed bar implies evaluation at b. The error in Eqs. 
(5.2.9) and (5.2.10) arises from: (1) the truncation of higher order 

moments and (2) the discretization of the random field b(x) by the 
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finite vector b. If the randomness in b(x) is small, then the first 
error will be small for a smooth function and the second-moment analyis 
is applicable. The error due to discretization in Eqs. (5.2.6) and 
(5.2.8) can first be studied to provide an insight of the discretization 
accuracy . 

The discretization error of the covariance field is defined by the 
L2 norm 

E 2 - / [Cov E (b(x k ),b(x A )) - CoVpCbCx^) ,b(x^))]^dft (5.2.12) 

where Covg and Cov^ represent the exact and discretized covariances. 

The exact covariance is calculated from the given function for the mean 
E[b], coefficient of variation u and the autocorrelation R as follows: 

Cov E (b(x k ) ,b(x 4 )) ■ [Var(b(x k ))Var(b(x 4 ))] 1 ^ 2 R(b(x k ),b(x fc )) (5.2.13) 

where 


Var(b(x k )) - (a(b(x k ))E[b(x k )] ) 2 (5.2.14) 


The discretized covariance between any two points x, and x is obtained 

~k ~2, 

from: 




lj-1 


(5.2.15) 
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and 


Cov(b i ,b j ) - [Var(b 1 )Var(b j )] 1/2 R(b i ,b j ) (5.2.16) 

where b. are the discretized points of b(x), corresponding to x , i - 1, 
..., q. For a beam with a random field along the x-axis (Fig. lb), the 
logarithmic plot of E against q is given in Fig. 2a, and the rate of 
convergence is found to be 1.325 for nearly all q between 4 and 64. 

When the random field discretization is coupled with an FEM 
discretization, as in PFEM [13], q need not be equal to the number of 
finite elements NUMEL and the shape functions N. (x) need not be the same 
as the finite element interpolants for the displacement field. 
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5.3 Transformation Procedure for Computational Efficiency 

The mean and covariance can be obtained from Eqs. (5.2.11) and 
(5.2.12). However, the number of derivatives to be evaluated Is 
proportional to q(q+l)/2. This arises from the double summations in i 
and j. To reduce the number of computations, the full covariance matrix 
Cov(b^,bj) is transformed to a diagonal variance matrix Var(c^,Cj) such 
that 


Var(c i ,Cj) - 0 for i j 4 j 


(5.3.1) 


and 


Var(c^,Cj) - Var(c^) for i ■ j 


(5.3.2) 


Therefore, the number of derivative evaluations is proportional to q. 
The above is achieved through the elgenproblem: 


Gfp * i|/A (5.3.3) 

where the G and A matrices denote Cov(b i ,bj) and Var(c^,Cj), 
respectively; ip is a constant q x q fundamental matrix with the 
following properties: 



(5.3.4) 
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k “ t Sfc 


(5.3.5) 


and 


b 


jj»c or 



(5.3.6) 


1 is Che q x q Identity matrix and c is the transformed q x 1 vector of 
random variables. 

With Eqs. (5.3*5) and (5.3.6), the mixed derivatives appearing in 
Eq. (5.2.11) reduce to second derivatives and Cov(b^,bj) reduces to 
Var(c^) : 


E[S] 



q a 2 ? 

I ■*— !• Var(c. ) 

i-1 3c^ 


(5.3.7) 


and 


q 3 S S S 

c5 - 3 - 8> 

Thus, the discretized random vector b is transformed to an uncorrelated 
random vector c, with the variance of c as the eigenvalues of G in Eq. 
(5.3.3). 

In the numerical examples, an exponentially decaying 
autocorrelation in one-dimension is assumed with various correlation 
lengths *X' (i.e., the length at which the autocorrelation drops to 
0.37, see Fig. lb). It is observed that for one-dimensional random 
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fields, as X increases from zero to a large value, the number of largest 
eigenvalues N, N < q necessary to evaluate the mean and covariance in 
Gqs. (5.3.7) and (5.3.8) to a specified accuracy, decreases from q to 
1. When X is zero the random field is uncorrelated and all q 
eigenvalues are dominant. When the field is uncorrelated, all q random 
variables are necessary to represent the randomness of the field. 

As X increases the number of dominant eigenvalues decreases. 

Eventually, for a very large X the random field is closely correlated 
and there is just one dominant eigenvalue. When the field is closely 
correlated, only one random variable, corresponding to the largest 
eigenvalue, is sufficient to represent the randomness of the field. 

This feature, when present, can easily be exploited to reduce the 
number of computations. The value of N can be chosen based on the 
distribution of the eigenvalues before solving the PFEM equations. The 
eigenvalues here can be Interpreted as weighting factors for the 
corresponding mode shapes necessary to represent the covariance 
structure; a large eigenvalue means a dominant mode and vice versa. The 
eigenvalue distribution and the mode shapes are depicted in Figs. 2b, 3a 
and 3b, for a numerical example. Results of the eigenvalue distribution 
and selection of N, for a beam problem and a bar problem, are discussed 


in the next section 
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5.4 Results and Discussion 

(1) Elastic/Plastic Beam with Yield Stress as a Random Field 

The problem statement Is depicted in Fig. lb. The yield stress is 
assumed to be a function of the position along the length of the beam 
only. The Gaussian random field , which is the yield stress, is 
discretized so that q - 16 (NUMEL * 64). 4-node continuum elements are 
employed. The coefficient of variation of the yield stress is assumed 
to be 0.10. The static response, as a function of the loading, is 
calculated by an implicit algorithm. 

The mean displacement, the variance of the displacement, the meah 
bending stress and the variance of the bending stress are shown in Figs. 
4a to 4d. The coefficient of variation of the displacement at the free 
end is found to be 0.069 and that of the stress at the fixed end is 
0.087. The results are compared with those of a Monte Carlo simulation 
with 400 realizations and they are in excellent agreement. 

The convergence of the random field discretization error, as 
defined in Eq. (5.2.13), is plotted in Fig. 2a. The rate of convergence 
is found to be 1.325. The eigenvalues of the covariance matrix are 
plotted in Fig. 2b. Based on the distribution, 4 out of the 16 largest 
eigenvalues were chosen. The mode-shapes, corresponding to these 4 
largest and 4 smallest eigenvalues, are shown in Figs. 3a and 3b, 
respectively. The latter, clearly, play no role in representing a 
smooth autocorrelation, as is assumed here; if the field is highly 
uncorrelated these modes will be necessary. This resulted in a 95% 
accuracy in the variance of the stress at the wall (Fig. 2d). 
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The exacC autocorrelation and the discretized autocorrelation for 
the Monte Carlo simulation of 400 realizations are compared in Fig. 2c; 
the autocorrelation is along the length of the beam, w.r.t. the yield 
stress at the wall. This amply demonstrates that this sample size would 
be sufficient to bring out the response correlation characteristics, and 
that the first-order filter captures the correlation characteristics 
quite well. 

The spatial autocorrelation of the displacements at two different 
loads, along the length of the beam, are depicted in Figs. 5a and 5b. 

The spatial autocorrelation of the stresses along the length of the 
beam, at these loads, are depicted in Figs. 5c and 5d. The displacement 
autocorrelation is w.r.t. the free end displacement and the stress 
autocorrelation is w.r.t. the wall stress. In the first loading, 4 
element layers (16 elements) near the fixed end are yielded and in the 
second loading 10 element layers (64 elements) are yielded. The 
displacements along the length of the beam show almost complete 
correlation with one another, irrespective of the correlation 
characteristics of the yield stress. The stresses, because of their 
direct dependence on the material properties, exhibit a varying 
autocorrelation along the length of the beam just like the random yield 
stress. Interestingly, the results of stress autocorrelation by PFEM 
and MCS are smooth and in good agreement in those elements that have 
yielded (in the mean sense) (Figs. 5c and 5d). In other elements the 
results show disagreement and oscillations in both PFEM and MCS 
results. In many of these elements, the mean stress is well below the 
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mean yield stress and so the randomness of the yield stress has a very 
insignificant effect on the stress values, as measured by the stress 
variances. Figs, 5c and 5d show that the stresses in the unyielded 
portion of the beam are highly uncorrelated to those in the yielded 
position. It is interesting to note that although the material law is 
highly nonlinear, the second order moment method which underlies PFEM 
agrees very well with Monte Carlo simulations, 

(2) Elastic/Plastic Bar with Plastic Modulus as a Random Field 

The problem statement is depicted in Fig. 6a, The plastic modulus 
&P is assumed to be a Gaussian random field along the length of the 
bar. The material is assumed to be elastic-plastic with isotropic 
hardening. As can be seen from Fig, 6a, the yield stress is spatially a 
linear function for the mean and an exponential function for the 
autocorrelation. The coefficient of variation is assumed to be 0,10 and 
the random field is discretized so that q - NUMEL * 32. The 
probabilistic equations are solved by the explicit predictor algorithm 
[12] with a slight numerical damping (y * 0,55). A near-critical time 
step (At * 0.000455) is used to keep the number of time steps minimal, 
subject to the stability conditions. 

The mean and the variance of the displacement at the free end are 
shown in Figs. 6b and 6c. The coefficient of variation of the 
displacement at the free end is found to be ~0,05, The results are 
compared with those of a Monte Carlo simulation with 400 realizations 
and they are in excellent agreement. For both examples, the PFEM needed 
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much less computer time than the MCS. Also, only the largest 8 of the 
32 eigenvalues are found to be sufficient to predict the displacement 
mean and covariance with a 99 % accuracy. 

In second moment PFEM, the superposition of the covariances of the 
response for two different, uncorrelated (to each other) random fields 
in a structure is the same as when both the random fields are present 
simultaneously. For example, the results of a bar with only the yield 
stress as the random field and those of the bar with only the plastic 
modulus as the random field can be superposed. The summed results will 
be the same as that of the bar with both the random fields present. 

When ’N' random fields are present, they are divided into 'n' groups (n 
< N) such that the fields within a group are correlated to one another 
and uncorrelated to those in the other groups. The PFEM results for 
each of these groups can then be simply summed, as in the case of 2 
uncorrelated random fields. This is, of course, not possible in 
simulation where the entire calculations have to be repeated. For the 
purpose of probabilistic analysis in multi-component systems, this is an 
added advantage of PFEM over the simulation methods. 
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Some Negligible Modes 16 Most Dominant Modes 
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Comparison of Mean Displacement 
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Comparison of Variance of Displacement 





CHAPTER 6 


SUMMARY AND CONCLUSIONS 

The theme of this study was to develop and apply efficient 
probabilistic finite element methods for various classes of problems in 
structural and solid mechanics. In nonlinear problems the main issue 
addressed was the evaluation of the higher order derivatives of the 
stresses and the internal nodal forces with respect to the random 
variables (Chapter 2). It was shown that finite-differencing was a 
fairly accurate way of approximating these derivatives. Applications in 
truss structures were studied and the results agree with those of Monte 
Carlo simulation and Hermit e-Gauss ian quadrature. It was also 
discovered that the higher-order equations involve secular terms which 
caused the solutions to deteriorate with time. Damping did not mitigate 
the problem and other remedies were suggested to eliminate secular 
terms. Secular terms, however, do not arise when only the external 
forces are random. The spatial discretization procedure of the random 
field resulting in the mean vector and the covariance matrix was 
outlined in Chapter 3, along with a simple check on the discretization 
accuracy. The popular eigenvalue transformation technique to obtain a 
vector of uncorrelated variables has also been implemented in PFEM. 

This reduced the computations from a quadratic to a linear dependence on 
the number of random variables. Furthermore, it was observed that a 
reduced set of the uncorrelated variables was sufficient to model the 
randomness. Each variable represented a mode of correlation and 
depending on the strength of correlation adequate number of modes had to 
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be Included. In general, the higher the correlation the lesser was the 
number of modes that were required. The eigenvalues of the 
transformation were weighting factors In the superposition of modes. 
Whereas in Chapters 2 and 3 the PFEM equations were derived from the 
equations of motion, it was demonstrated in Chapter 4 that they can be 
quite as easily derived from variational principles. For linear 
contlnua, the equations were derived from the potential energy 
variational principle and for nonlinear contlnua undergoing large 
deformation the corresponding equations were derived from the principle 
of virtual work with appropriate stress and strain measures. Also, for 
elastoplastic materials, a direct method of evaluating the stress 
derivatives was outlined* An important advantage in deriving the PFEM 
equations from variational principles is the ease of incorporating the 
random geometry. 

Various applications were studied in truss structures, bar, beams 
and plates and the results compare favorably with those of Monte Carlo 
simulation. It was noticed that the second order terms have negligible 
contribution to the mean response. It was also noticed that the 
strength of correlation did not affect the magnitude of the response 
moments appreciably. However, this may not be the case when the 
structural reliability has to be calculated. With the use of adjoint 
method it was demonstrated that the response moments can be selectively 
computed in a desired portion of the structure. Some applications in 
elastic-plastic dynamics were studied in detail in Chapter 5. Most 
uncertainties seem to linearly affect the response (i.e., a 10% c.o.v. 
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of a random property gives to not more than 10Z c.o.v. of the 
displacement or stress). It was also noticed the displacements in any 
structure were always perfectly correlated regardless of the correlation 
level among the random variables. This was not the case for stresses, 
however. Unlike stress, which can exhibit sudden variation in a 
structure as in stress concentration, the displacement field in a 
structure is a smooth function by nature. Alternatively, from the 
finite element context, the external force is obtained by integrating 
the stresses over the domain and then the displacements are obtained. 
This integration process results in a smooth displacement field. 

Based on the work in this research, the following suggestions are 
made for further work: 

1. Methodologies for modeling mixed random boundary conditions, with 

PFEM. Direchlet and Neumann boundary conditions, with random 
magnitudes, can be easily incorporated in PFEM as random external 
forces. However, mixed boundary conditions are not so easily 
modeled. Example: A beam may not be fully clamped or fixed. 

2. Second— moment based reliability techniques are available for mostly 
linear systems. The major issue in nonlinear systems is the 
computation of the response gradients. By incorporating some of the 
techniques outlined in this report, this computation may be done, 
accurately and efficiently. In limited situations, the response 
moments themselves may be used to calculate reliability. 

3. As PFEM involves solution of a set of independent first-order 
equations, in addition to the zeroth and the second-order equations. 
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Che flrsC order equations can be processed in parallel. The 
solution sequence would be: solution of the zeroth-order equation 

first, parallel solution of the first-order equations next and 
finally the solution of a single second-order equation. In solving 
the higher order equations, the already decomposed stiffness matrix 


can be made use of 
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APPENDIX A: REVIEW OF MONTE CARLO SIMULATION METHODS 

AND HERMITE-GAUSS QUADRATURE SCHEMES 

A.l Monte Carlo Simulation Methods 

Various Monte Carlo Simulation techniques are now available. In 
our analysis, the "simple” Monte Carlo method is used. The important 
feature about the Monte Carlo method is its flexibility. In other 
words, the computational procedures are the same irrespective of whether 
the model is linear or nonlinear so long as the solution can be obtained 
from the governing equation. 

The main idea behind Monte Carlo simulation methods is to randomly 

generate values of the random variables subject to the probability 

density function and to calculate the output corresponding to these 

values. From this set of output, the probabilistic distribution 

properties, such as the mean and variance, are statistically estimated. 

In the analysis of the two degrees of freedom probabilistic system, 

a normal random generator is used. This normal random number generator, 

RANF, is available on the Northwestern University CDC system. It has 

been well tested and for large sample size, the distribution is close to 

normal. This "closeness” can also be estimated by the so called Central 

Limit Theorem, which states as follows: "if a population has a finite 
2 

variance a and mean value y , then the distribution of the sample mean 

2 

approaches the normal distribution with variance a /n and mean 4 as the 
sample size n increases." The sample size used in this analysis is 
400. Both the spring constants and are randomly generated and the 


corresponding displacement solutions are calculated using the exact 
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deterministic solutions. 

A. 2 Implicit Time Integration with Hermlte Gauss Quadrature Scheme 
Let us consider a linear two degrees of freedom probabilistic 
system 


M a + K d 
~ ~n+l ~ ~n+l 


-n+1 


where 


(A.2.1) 



(A.2.2) 

n Is the time step number and the initial conditions d , v and a are 

~o ~o ~o 

given. At each time step n, Gq. (A.l) is solved by the Newmark-3 

algorithm with 0 * 0.25, y ■ 0.5 and At * 0.02 T . where T . is the 

1 min min 

smallest fundamental period. The finite difference matrix equation can 
be shown to be 

K* ff d . - F* ff (A.2.3) 

~ ~n+ i ~ 


where 
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(A.2.4) 


K* ff - M + SAt 2 K 

/w /%» ~ 


(A.2.5) 
(A. 2.6) 

and 


P* ff =■ 0At 2 F + M 3 
- M ~n+l ~ ~n+l 


d - d + At V + (~ - s) At 2 a 
~n+l ~n ~n Z 


V, - v + At (1 - y) a 
~a+l ~n ~n 

Once d Q+ j is determined from Eq. (A. 2. 3), i.e., 

V. - (K *”)' 1 P' ff 

a n+ j and v Q+ ^ can be determined as follows 

2 , 1+1 ’ < 2 , 1+1 " 3 n + l > /61t2 

and 


(A.2.7) 


(A. 2.8) 


(A. 2.9) 


v 

~n+l 


•Vi * 


yAt a 


n+1 


(A.2.10) 


The solution procedures are then repeated with n replaced by n + 1 
until nAt greater or equal to a desired time. Since and ^ are 
random variables, therefore, by Eq. (A.8) is ’’implicitly" a 
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function of and K^> Using the basic definitions of mean and 

variance, the expected value of d , , is 

~n+i 


<A»ll ' i i 2a.l<Y V \ (K 2> dK l dK 2 (A - 2 ' U > 


where f (K ) and f (K ) are the probability density functions for 

aV^ X 1^2 4 

and respectively. In writing Eq. (A.2.11), the assumption that 

is uncorrelated to has been employed. Once the expected value is 

evaluated, the variance of d can then be computed according to 

~n+l 


v * r (ia.l) - 'tS.-ll - ( E l d o.l ]) 2 


(A.2.12) 


The Hermite Gauss Quadrature scheme is to approximate the double 
integrals which appears in Eq. (A.2.11) by 


nl n2 


^n+J * ^ W i W j ~n+l (K l’ K 2* f K 1 (K l ) f K^ K 2^ 


(A.2.13) 


where nl and n2 are the number of integration points for and 
respectively, w^ and w are their corresponding weights. As might be 
figured from above equation, if the number of random variables is m, the 
number of simulations N is 


N«nl*n2*....*nm 


(A.2.14) 


and N grows exponentially. Therefore, unless the number of random 
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variables Is small, this method Is not recommended. However, If a 
physical situation dictates that some of the random variables can be 
excluded In the calculations, N can then be reduced significantly. 

Under this circumstance, the Hermite Gauss Quadrature Scheme can be an 
efficient and accurate method. An example of this practical situation 
has been demonstrated in Chapter 2, Section 4. In the analysis of the 
ten-bar structure, It was predetermined that only four of the ten bars 
will yield and the yield stresses are chosen as normal random variables. 

If 3 points are used in evaluating each normally distributed in 
Eq. (A. 2. 13), the weight and quadrature points are 

w ± - (1/6, 4/6, 1/6) (A.2. 15a) 


and 


- (u - 3cr , u , u + 3<j) (A.2. 15b) 

respectively. And for this "bias" integration procedure, the number of 
simulations required (for each time step n) is 

N-3*3*3*3-81 (A. 2. 16a) 

Whereas, if one cannot observe apriori that six of the ten bars will not 
yield, the number of simulations required becomes 


N - 3 10 - 59049 
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(A. 2. 16b) 


which makes this a handicapped method 


APPENDIX B: RESONANT EXCITATION OF RESPONSE SENSITIVITIES 


The equation of motion and the sensitivity equation for a single- 
degree— of-freedom spring-mass-damper system are 


Mx + Cx + Kx- F(t) 


(B.l) 




(B.2) 


where 


F(t) 


3F 3M *• 3C . 3K 
3b ” 3b x " 3b x “ 3b 


(B.3) 


M, C, K are assumed to be dependent on the parameter b; we are 
interested in the sensitivity of the response x(t) to this parameter 
b. Let F(t) be such that x(t) is stable. Under this condition it is 
shown below that the response sensitivity ^ is resonantly excited. 

The damped natural frequencies of the system Eq. (B.l) and the 

A 

sensitivity Eq. (B.2) are the same. The excitation F(t) involves 

x, x and x in Eq. (B.2) and is, therefore, a resonant excitation. Thus 

approximations for x(t), such as 

2 

x(b + Ab, t) * x (b ft> t) + (|f) Ab + i (^-f) Ab 2 (B.4) 

d b-b 3b b-b 

o o 
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at b for any small Interval Ab, are valid only for a short duration 
o 

and the accuracy deteriorates rapidly thereafter. 

Since the PFEM equations (2.7) in Chapter 2 use the first and 
second-order response sensitivities, they are valid for a short duration 
only. A similar phenomenon is also observed in the transient response 
of nonlinear structures. A possible explanation for this phenomenon is 
that the time 't* has a multiplying effect on the interval 'Ab' in the 
second and third terms in Eq. (B.4) and this results in the 
deteriorating accuracy. However, the PFEM equations developed in 
Section 2 of Chapter 2 are suitable for application when one is 
interested in a short-time history e.g. , the response due to an 
impulsive load. 



APPENDIX C: NUMERICAL ALGORITHM FOR STRESS 

AND DISPLACEMENT DERIVATIVES 


For materials with random elastoplastic properties, such as yield 
stress or plastic modulus, the finite-difference based derivatives (Eq. 
(4.4.15b), Chapter 4) can be evaluated explicitly without recourse to 
the solution of the equilibrium equation at each finite-differencing 
point; the stiffness matrix, corresponding to the mean configuration of 
the random properties, needs to be decomposed only once* This 
configuration is represented by the vector c of size q. Subsequently, 
the displacement derivatives can be calculated by forward reduction and 
back substitution in Eq. (4.3.15a), Chapter 4* These computations are 
done in conjunction with the radial return method proposed in Ref. [1] 
and as implemented in Ref. [2]. 

For each element integration point, for every equilibrium iteration 
'v f in a given load step ’m 1 , the following quantities are stored: 

7 * : stress 

r J center of the von Mises yield surface 
k s radius of the yield surface 
C T : tangent constitutive matrix 
where the superposed ” implies quantities evaluated at c. 

FLOWCHART 

Part I Radial return for mean stresses and dl placements 
Begin loop on load steps, m 

Begin loop on iterations, v 
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Begin loop on elements 


Begin loop on integration points 


a) Pick up c v , o'* , k'* , C 

b) Compute the elastic trial stress Increment aT'* and update the stress 

— v+1 


~ V , — V 

t + At 


to obtain the elastic trial stress , : x . , 

~trial ~trial 

c) Compute: £+J al - deviatoric part of (7£} al - a V ) 

d) Compute: A - 3J_ - (k 

z n 

where Is the second invariant of stress and $ is the yield 
function. 

e) If J < 0, 7 V+1 - T^ al ; go to step h. 

If > 0, compute the plastic strain Increment Ad 1 *: 

Ad 9 - (>/3^ - 1T , )/(3G + B) 

where G and B are the shear modulus and plastic modulus. 


(c.l) 

(C.2) 


(C.3) 


(C.4) 


respectively. 

B .Jfl 

B E-E_ 


' B* is further defined as: 


(C. 5) 


where E and E^ are the elastic and tangent moduli, respectively. 


f) compute: 


AS - 3G Ad 

(C.6) 

Ao - (1 - B)B AdP/V^ 

(C.7) 

Correct stresses by radial return: 


—v+1 —v+1 —v+1 

~ " Xtrial 8 2trial 

(C.8) 

—v+1 —v , . —v+1 

2 * 2 + A* Atrial 

(C.9) 

k v+1 - i?* + BB Ad 9 

(C.10) 
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h) Compute the tangent constitutive matrix C_ and assemble the tangent 

~T 

stiffness matrix 

End loop on Integration points 

End loop on elements 

End loop on iterations until equilibrium is satisfied. The 

stresses, strains and displacements, T , , T , and d" , are 

~nrt-l miH- 1 ~m+l 

obtained from the final iteration. 


Part II Displacement derivatives 

The following quantities are also stored for each element 
integration point, in a given load step m, for each 'i' representing a 
random variable c^: 


+i -i - 

t , x : stresses, near the mean c 
~m ~m ** 


a + *, a * : yield surface centers near the mean c 
~m ~m +* 


k +i , k ^ : yield surface radii, near the mean c. 
mm 


C+ 1 , : tangent constitutive matrices, near the mean c, 


where 


+i — i T 

c * c + (0, ... , Ac 0) 
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(C.ll) 


c * c “ (0 j * • * > Ac ) • . . > 0) 


(C. 12) 


and 


Ac 1, - Be 1 , 0.01 < 0 < 0.05 


(C.13) 


Begin loop on for each random variable 


Begin loop on elements 


Begin loop on Integration points 


1) Compute the strain Increment Ae as 

~{Q 

* T , t - 7 

~nH* i <»»m 


(C.14) 


... _. . +1 +1+1 . +1. 

11) Pick up t . a . k . (C ) 

~ra ~m m ^ i m 

ill) Compute the elastic trial stress (x + ^,) tr ^ a1 ’ 

~nH-l 

lv) Repeat steps (c) through (h) to finally obtain 


+i 


+i 


r n k"^ 1 fC +1 l 

~mfl ’ 2m+r m+1’ ~T ; m+l 


v) Pick up t * . o * , k 1 , (C_) 1 and repeat steps (iii) and (iv) to 
~m ~m m T m 

.Iso obtain Cc; i ) ntfl 
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vi) Compute and store the first order finite-difference derivatives: 


3t . , 
~m+l 

3c. 


+i 


-i 


~m+ 1 ~m+ 1 

2Ac. 


(C. 15) 


3C C? - C" 1 

and - — a — sr 

3c i 2Ac^ 


(C. 16) 


vii) Compute and assemble the forcing term f^ in Eq. (4.3.15b), 
Chapter 4; note that 


4 


• 3 Im+l 


3c. 


(C. 17) 


in the equation. 

End loop on integration points 


End loop on elements 

viii) Solve Eq. (4.3.15a), Chapter 4 to obtain the displacement 

3d., 

/v*QTT I , 

derivatives and store. 

3c, 


Part III Stress derivatives 


Begin loop on elements 

Begin loop on integration points 

ix) Compute stress derivatives as: 


[7 i . !~2±I + c b r-^2i±') 

L ~m+l J Ci 3c £ ~T~ l 3c t } 


(C. 18) 
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End loop on integration points 


End loop on elements 


End loop on random variables 


End loop on load steps 



